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Investigation of Lateral Dynamic Stability 
in the XB-47 Airplane 


ROLAND J. WHITE* 


Boeing Airplane Company 


ABSTRACT 


This paper presents the method and analysis used to improve 
the lateral stability of the XB-47, which was a part of the U.S. 
Air Force development program on this airplane. Early flights 
indicated a tendency for the airplane to ‘“‘Dutch Roll’’ at the 
lower indicated air speeds. Elimination of this flight condition 
has been accomplished by use of a servomechanism termed a 
“yaw damper,” consisting of a yaw rate gyro pickup that oper- 
ates the rudder to improve the airplane damping in yaw. The 
resulting flying qualities are believed to be considerably better 
than would be possible with other means of improving the dy- 
Problems encountered during the design and 
A method of analysis based upon the 


namic stability. 
flight tests are discussed. 
airplane inverse frequency response characteristics has been used, 
from which the airplane stability has been obtained. This 
method has been extended to include the effect of airplane body 
side bending which influences the dynamic stability at higher 


indicated air speeds. 
INTRODUCTION 


_— LATERAL STABILITY of the XB-47 was investi- 
gated in the early design stage, using stability 
equations including the effect of the airplane product 
of inertia.! The analysis was based upon static sta- 
bility coefficients obtained from wind-tunnel tests, 
along with calculated damping coefficients. 

Early flight tests indicated that the calculations were 
generally correct but that the damping in yaw, par- 
ticularly at low speeds and at altitude, was less than 
the pilots desired. The most important factor was 
that low damping could not be tolerated even though 
the periods were long. For example, during the land- 
ing approach the pilot’s attention was required to pre- 
vent a Dutch Roll from occurring, which would be a 
handicap during blind landing approaches. The prin- 
‘Sipal objection of the pilots to the poor damping 
stemmed from the rolling part of the motion. This was 


Presented at the Design Session, Annual Summer Meeting, 
LAS., Los Angeles, July 21-22, 1949. 
* Aerodynamicist. 


later confirmed by pilot checks on an airplane having 
poor directional damping but low dihedral effect, 
causing only a small rolling. While the rolling part of 
the motion is more apparent to the pilot, it should be 
noted that, even if rolling were not present, the yawing 
motion would still be poorly damped. This condition 
is characteristic of most all high-performance airplanes 
having high wing loadings, which results in a high 
“Relative Density Factor” » = (W/g)/pSb, a term used 
in airplane stability analysis to express the ratio of the 
airplane mass to the mass of air enclosed in a volume 
equal to the wing area times a characteristic length 
of the airplane. It can be shown that the damping ratio 
is proportional to V ‘) /y when the airplane is considered 
as a simple vibrating system. Therefore, the number 
of cycles to damp to half amplitude is a function of uy. 
Since yw increases with altitude and wing loading, it 
should be expected that improvements to power plants 
and aerodynamics, permitting higher loadings and alti- 
tudes, should cause a decrease in damping, assuming 
that the general proportions of the airplane are not 
changed. Fig. 1 shows the relative density of a num- 
ber of airplanes at their service ceiling, indicating the 
increase in wu that has taken place in recent years. Some 
means must be found to improve the damping on all 
future high-performance airplanes, since the damping 
is already marginal in many cases. 

Following the early flight tests, during which a more 
clear understanding of the desired flying qualities was 
obtained, the stability was further analyzed to see 
which feature contributed to poor damping and to ex- 
amine means of altering the general configuration to 
improve the damping to an acceptable condition. The 
dihedral effect resulting from the wing sweepback was 
first studied and found to reduce the damping some- 
what, although the total magnitude of this effect was 
small compared to the desired increase in damping. 
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Fic. 1. Effect of time on relative density factor. 
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stability. 


It was also realized that if sufficient negative di- 
hedral could have been obtained for the flaps-down 
approach condition, this would have given marginal 
static sideslip stability in the flaps-up condition at high 
speeds. Fig. 2 shows the effect of reducing the wing 
dihedral during the approach condition. Here, a large 
change was out of the question because of ground clear- 
ance. This same figure shows the effect of increasing 
the vertical tail area, indicating that a vertical tail of 
twice the present size would be necessary if damping to 
half amplitude in 1 cycle was desired. 

While increasing the tail size increases the damping, 
it is problematical if any gain is made in the smooth- 
flying qualities when the vertical tail is increased above 
a certain size. This is because the excitation of the 


airplane due to side gusts will increase proportional ty 
the increase in vertical tail size; hence, while the damp. 
ing may be greater, so will the disturbances, Tesulting 
in a doubtful improvement in the airplane riding cop. 
fort. It is probable that the usual tail volume q. 
efficients that satisfy the requirements for unsymmet. 
rical power, static stability in sideslip, improved stal}. 
ing qualities, etc., should determine the tail size, ang 
some other means should be used if large increases jy 
the damping are required to obtain the desired riding 
comfort. 

Inasmuch as a rate gyro servomechanism was most 
easily assembled, consideration was given to improving 
the damping in the yaw axis. 

The improvement of airplane stability by automati- 
sally moving a control surface in response to an airplane 
disturbance has been used in many instances. These 
devices range from springs, bob weights, and aerody- 
namic vanes to the present yaw damper or similar 
servomechanism devices—the latter being a more suit- 
able, and a positive means, when power operated irte- 
versible controls are used. The distinction between 
these devices and the normal autopilot is that they im- 
prove the airplane flying qualities so that the human 
pilot may more easily fly the airplane, whereas the 
normal autopilot adds direction sensing and has com- 
plete authority of the airplane flight path. 

The use of a rate gyro to improve the damping in 
yaw was decided upon from practical installation con- 
siderations. This was initially thought to be a new 
idea, but later it was realized that this concept has 
been used at least in one known case.? The analysis 
in this paper will be limited to this type of control. 

The method of analysis used has been useful in de- 
veloping the yaw damper. Unfortunately, in several 
cases the analysis lagged the flight experience. How- 
ever, the analysis has served to give a clear picture of 
the damper operation and was invaluable in making 
laboratory and flight-test changes. 


DESCRIPTION OF YAW DAMPER 


A general view of the XB-47 airplane is shown i 
Fig. 3. All control surfaces are power-operated and 
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Fic. 3. General arrangement of XB-47. 
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Fic. 4. Schematic diagram of rudder system with yaw damper. 


are irreversible for their normal operation with arti- 
ficial pilot “feel.” A sufficient amount of aerodynamic 
balance is provided for each control surface, so that 
control forces with the servo-off, while higher than 
desired for normal flying, are reasonable for emergency 
operation. As this paper is concerned with the rudder 
control system, a detailed description of this system is 
given. 

A schematic diagram of the rudder controls located 
at the tail of the airplane is shown in Fig. 4. Their op- 
eration will be described first without consideration of 
the directional control damper. During ‘‘boost-on” 
operation, movement of the rudder pedal produces a 
rotation of the rudder torque tube, which is resisted by 
a torque at its lower end, developed by the artificial 
feel spring. The control arm at the upper end of the 
torque tube moves the push-pull tube to position the 
hydraulic boost control valve, which requires only a 
friction-resisting force of 3 to 5 Ibs. In turn, the hy- 
draulic boost cylinder moves the rudder until the follow- 
up link centers the boost valve at the desired rudder 
position. Pilot pedal forces are resisted at the lower 
end of the vertical torque tube by a force gradient, 
which is essentially linear with rudder deflection and 
proportional to the dynamic pressure. This is ac- 
Point A of the 
torque-tube trim yolk remains fixed and movement of 


complished in the following manner: 


the torque tube rotates point B, causing point C als» 
to rotate. 
D. The artificial feel spring consists of a piston and 


This produces a sideward movement of point 


fabric seal, providing a tension force in the cable at- 
tached at point D, which is located at the end of an 
arm that is free to rotate about the torque tube. Side- 
ward movement of point D provides a force proportional 
and to the deflection of point D. When the 
pilot pedal force is released, point D centers itself and 
When it is desired 


to a 


point C becomes essentially fixed. 


to trim the rudder with ‘“‘boost-on,’’ movement of the 
trim control operates the jack screw to move point A. 
Here, the only resisting torque is that of the pedal con- 
trol cables and the boost valve; hence, point B will 


move because the feel spring force keeps point C fixed. 
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In this manner the boost valve is positioned, causing 
a cliange in the rudder position. During this opera- 
tion parallel cables to the rudder trim tab (not shown) 
will tend to trim the rudder aerodynamically, so that, 
should the boost be turned off, the rudder will still 
remain closely in trim. When the boost is turned off, 
loss of hydraulic pressure relieves the air pressure in the 
“q’’ spring container and locks the boost follow-up 
motion so that the push-pull tube operates the rudder 
directly with no lost motion. 

The directional control yaw damper has a rate gyro 
designed to develop an electrical signal voltage pro- 
portional to the rate of yaw of the airplane. | This volt- 
age then passes through a servo-amplifier to eperate an 
electric push rod positioning motor. The motor oper- 
ates the push-pull tube shown in Fig. 4 to extend or 
shorten it by an amount proportional to the electrical 
signal from the rate gyro. The push-pull tube at the 
torque-tube end is restrained from moving by the 
“q’’ spring; hence, the change in length positions the 
boost valve, causing the rudder to move. This push 
rod only operates the boost hydraulic valve; hence, 
the operation loads are small. The full extension and 
shortening of the push-pull tube is +0.40 in. and will 
produce a +5° of rudder movement. This has been 
found sufficient for the maximum oscillation amplitudes 
that the pilot would expect the airplane to damp with 
free controls. The operation of the rudder occurs with- 
out motion of the rudder pedals; hence, the pilot is not 
aware of the yaw damper operation, except in noting 
the more desirable flying qualities of the airplane. 
Since the yaw damper adds a differential signal to the 





Rudder boost push rod. 


Fic. 5. 
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rudder, the pilot can overpower the yaw damper. An 
angular yaw velocity of the airplane will cause the gyro 
to precess about the transverse pivot axis, which is 
parallel to the lateral axis of the airplane, until the 
movement is balanced by springs restraining the mo- 
tion. The deflection of the gyro gymbal develops a 
signal voltage in the yaw damper servo system, which 
acts to position the rudder. 

The original yaw damper servo system was made 
up, for a brief time, out of surplus electrical parts. 
This was for the purpose of trying out the general idea, 
although it was anticipated that deficiencies in the 
response of the servo might occur. As a result of the 
flight and laboratory tests, various improvements were 
made to the system. Calculations in this paper are 
based on the yaw damper operating with the improved 
Fig. 5 shows a view of the rudder boost push 
The general 


system. 
rod and its installation in the vertical tail. 
operation of the rudder system may be visualized from 
the block diagram arrangement of Fig. 6, which shows 
how the yaw damper superimposes its motion with that 
of the pilot’s pedal to operate the rudder. 


RESULTS 


Flight-test graphic records of the airplane motion are 
shown in Figs. 7 and 8. Here the airplane is made to 
Dutch Roll with the damper-off showing the low damp- 
ing. The pilot then turns the damper on and the 
Dutch Roll rapidly disappears. Fig. 7 represents the 
approach flight condition and Fig. 8 is typical of high- 
altitude flight. 

The calculated periods and time to damp to half 
amplitude, along with several flight test check points, 
are shown in Fig. 9. These have been compared with 
the recent U.S. Air Force requirements, indicating that 
the airplane is unsatisfactory or marginal in several 
flight conditions without the yaw damper. Calcula- 
tions with the yaw damper operating show an effective 
improvement. Flight have indicated that 
higher damping is desirable; hence, in Fig. 9, a line 
marked ‘Desired by pilots’’ has been added. Fig. 10 
is a stability diagram in which the damping in yaw and 
This 


figure shows how the yaw damper increases C,,; to move 
d ty 


results 


airplane static yaw stability form the variables. 


the airplane point away from the oscillatory instability 
boundary, thereby improving both the Dutch Roll 
and spiral stability. 

The magnitude of roll to yaw amplitude ratio is of 
interest and is shown in Fig. 11. Here the roll to yaw 
amplitude ratio changes because of different rudder 
input frequencies, which, of course, equals the airplane 
frequency. At the resonant frequency wz of both the 
high- and low-speed cases shown, the amplitude ratio 
At higher frequencies, however, the roll 
For this reason, the effect of 


is near 1.9. 
effect becomes small. 
roll may be neglected where higher frequencies are con- 
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cerned, such as are encountered in certain types of 
hunting problems with automatic controls. 


FLIGHT-TEST EXPERIENCE 


The first low-speed flight test with the yaw damper 
was extremely successful. Here the yaw damper re. 
moved the Dutch Roll tendency, and no hunting of the 
airplane was experienced. This flight record is showy 
in Fig. 7. In the following flight, the yaw damper 
rudder stops were set at 59 = +3° instead of the pre. 
vious setting of 69 = The 3° 
to reduce the effectiveness of the yaw damper; hence, 


ee setting was found 
the stops were reset to 5° and have been held there 
since. 

One interesting aspect of the early flight tests was 
that a rheostat control was installed at the pilot's sta. 
tion to vary the sensitivity of the yaw damper. The 
sensitivity is defined as the amount of rudder movement 
produced by a steady unit rate of yaw. The pilot could 
vary this from 0.3 to 1.3 (degrees of rudder per degree 
per second yaw velocity). Invariably the pilot would 
report that the flying qualities were best when the 
sensitivity setting was maximum, giving maximum 
damping. It also appears that the pilot would like all 
airplane motions to appear ‘‘dead beat’’ to him. Actu- 
ally, this requires only about 20 per cent of critical 
damping. 

Initial flights were limited to an IAS = 200 mph 
until the damper operation was explored. When higher 
indicated air speeds were tried, a definite airplane 
hunting was experienced at a frequency greater than 
the natural airplane frequency. This was slightly 
noticeable and intermittent near IAS = 300 mph 
but became fully developed at IAS = 380 m.p.h. Cal- 
culations indicated that the servo system was not fast 
enough, causing the large phase losses to occur at too 
low a frequency, as will be shown later. The servo 
system was then improved with a faster motor that 
avoided the previous type of hunting. Another type 
of hunting then appeared at a frequency close to the aif- 
plane’s natural frequency. This hunting motion was 
suspected as being due to the phase lags developed at 
extremely low yaw amplitudes, which is a commot 
problem in the development of all servos. To improve 
this situation, the sensitivity of the internal electrical 
servo was increased, requiring additional complications 
to maintain stability of this part of the servo loop 
This made the yaw damper more sensitive at low ampl 
tudes and reduced the lag. Flight tests with this 
system showed that the system was improved to the 
point where it could be used satisfactorily in service 
Examination of the phase lags throughout the complete 
system, including the hydraulic surface boost, is neces 
sary before the system can be considered satisfactory 4 
the smaller amplitudes of motion. 

During the recent tests, a high-frequency hunting 
developed which was associated with body bending. 
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This was corrected by moving the rate gyro near the 
airplane center of gravity. A technical discussion of 
this will be given later. 

The pilots general acceptance of the yaw damper 
has been good. The fact that the yaw damper operates 
the rudder without transmitting its motion to the rudder 
pedals has been found to be favorable. During rapid 
rolls or entries into turns, the yaw damper has been 
beneficial. The only objectional features have been 
the airplane small hunting amplitudes and the increased 
pedal forces during the steady portion of the turns. 
Yaw damper modifications to eliminate the increased 
control forces during the turn will be discussed in a 


following section of this paper. 


METHOD OF DETERMINING AIRPLANE STABILITY FROM 
FREQUENCY RESPONSE CALCULATIONS 


The usual procedure for calculating airplane stability 
involves the solution of the airplane stability equations 
for the roots, which enables the period and time to 
damp to be evaluated. A similar procedure could be 
applied to the XB-47 airplane if an effective airplane 
damping coefficient due to the yaw damper were esti- 
mated. This method would not, however, conveniently 
predict any instability conditions that may exist for 
the yaw damper while in combination with the airplane; 
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hence, it is of little use for design purposes. For thc 
reason, the present analysis based upon a knowledge of 
the airplane and yaw damper steady-state frequency 
response characteristics has been used. By properly 
combining these frequency response curves, a method 
of estimating the period and time to damp of the ai. 
plane with or without the yaw damper is presented 
The basic principle of the analysis has been based upon 
use of the inverse frequency response technique.* 


(a) Airplane Frequency Response 


The airplane frequency response will first be dis. 
cussed by considering the airplane to have a single de. 
Forcing 
the airplane to yaw in a positive direction at low fre. 


gree of freedom in yaw, as shown in Fig. 12a. 


quencies requires the pilot to deflect the rudder in q 
negative direction according to the N.A.C.A. standard 
notation. 
the rudder will be considered as —6. 
expedient the N.A.C.A. standard notation can be used, 
and the resulting response curves will be recognizable 
If the rudder is con- 


For this reason a positive disturbance of 
By this simple 


to the servomechanism engineer. 
tinuously oscillated in a sinusoidal manner at a given 
frequency w and amplitude 69, the airplane will respond 
with a similar motion having the same frequency but, 


in ‘general, a different yaw amplitude Y. The airplane 
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Representation of airplane frequency response Ya. 
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motion will lag the rudder motion in that the maximum 
yaw amplitude will occur after the rudder has experi- 
enced its maximum amplitude. The relation of the 
yaw motion with respect to the rudder motion may be 
represented by a pair of vectors, as shown in Fig. 12b. 
When this vector combination is rotated at the angular 
velocity w for which they are drawn, their vertical com- 
ponents plotted against time provide a time-history 
plot of the airplane and rudder motion. Computing 
similar vectors for a series of frequencies will give a locus 
of the end points of the Yo/ — dy vectors forming a polar 
representation of the airplane frequency response. In 
Fig. 12c, this locus is calculated from the equations of 
motion of the airplane given by the airplane transfer 
function VY, considering a single degree of freedom. 
Equations for calculating Y, for three degrees of free- 
dom are given in Appendix A and were used in the 
actual stability calculations. The vector of maximum 
length of Fig. 12c corresponds to the resonant fre- 
quency of the airplane and is denoted as we. At this 
frequency the pilot obtains the greatest amount of yaw 
for the least amount of rudder motion. 

In the present analysis, it is more convenient to plot 
the airplane dynamic characteristic in terms of the 
amount of rudder amplitude and its phase relation 
necessary to sustain a continuous oscillation of the air- 
plane. This is called the “airplane inverse frequency 
response’ and is obtained by plotting the reciprocal 
of the length of each Y/—d9 vector using a negative 
value of the phase angle. Fig. 13 shows the Y, plot of 
Fig. 12c converted into the airplane inverse transfer 
function denoted as 1/Y,. Here the resonant fre- 
quency vector becomes the minimum length vector in 
the plot. For the single degree of freedom considered, 
this plot has a parabolic shape. The next problem is 
to consider the yaw damper frequency response and to 
see how it may modify the present airplane inverse 


frequency response. 


(b) Yaw Damper Frequency Response 


The airplane yawing motion causes the rate gyro to 
produce an electrical signal that feeds into an electro- 
servomechanism to produce a_ rudder 
movement. At low frequencies a positive airplane yaw 
directs the yaw damper to produce a positive rudder 
The rate gyro will basically pro- 


mechanical 


angle or a —(— 6p). 
duce a 90° phase lead and a rudder displacement signal 
This gives a fre- 
quency response curve VY; as shown in Fig. l4da. The 
resulting electrical signal then enters the electrome- 
chanical servo unit, which consists of an electro servo 
unit that operates the push-pull rod attached to the 
hydraulic boost, which in turn positions the rudder. 
For simplicity of this discussion, the frequency re- 
sponse of this combined unit is represented in Fig. 14b, 
where the transfer function Vy, is specified in a simple 
Here, the design variables are 


proportional to the yaw frequency. 


approximate manner. 
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the sensitivity G and the servo time constant 7's. The 
frequency response of the rudder motion due to a yaw 
signal is then the vector product of Y¢ VY and is obtained 
by multiplying the respective vector lengths and add- 
ing their phase angles to form the feedback frequency 
response Y, shown in Fig. l4c. This diagram repre- 
sents the combined yaw damper and boost frequency 
characteristics, and for the assumptions made it is a 
circle. 

The accurate estimation of VY, from the design char- 
acteristics of the servomechanism is usually difficult 
because of the complicated and nonlinear character- 
istics of such devices. Therefore, the airplane manu- 
facturer resorts to laboratory tests to measure the fre- 
quency response of existing automatic controls. This 
is done over the range of frequencies and amplitudes 
expected in flight. The results of typical laboratory 
tests will be shown later in this paper. 

(c) Combined Frequency Response of Airplane with 
Yaw Damper 

The yaw damper opposes any motion applied by the 

pilot to operate the rudder; hence, when the pilot ap- 
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plies a (—4) rudder deflection to give the airplane a 
positive yaw, the damper supplies the rudder with 
a —(—dépy) = dy motion. As the frequency response 

‘y of the yaw damper is the ratio of 5,/yo, it is possible 
to add each Y, vector to the 1/ ¥, vectors which repre- 
sent the input motion —60/p. The resulting polar 
plot then becomes the overall inverse frequency re- 
sponse of the airplane with the yaw damper. This is 
merely a vector superposition of the two motions oper- 
ating the rudder. The validity of this addition is given 
in Appendix C. Fig. 15 shows this vector addition, 
where the polar plot denoted as 1/ Yj is the airplane in- 
verse frequency response with the yaw damper oper- 
ating. 


(d) Conditions for Airplane Instability 


Suppose next that, upon adding the Y; vectors to the 
1/Y, vectors, the resulting 1/ Yj curve passed through 
the ‘origin. Then the reciprocal of the zero vector 
length at the frequency corresponding to the origin 
intersection would give an infinite airplane response in 
yaw, which is the condition that occurs for an undamped 
system when it is being excited at the resonant fre- 
quency. When the 1/¥q polar plot passes near the 
origin, the Y, vector is opposite and almost equal to 
the 1/ ¥,4 vector at some particular frequency. Under 
this condition the product of Y4 Ys approached 1.0 in 
magnitude, with the sum of the phase angles approach- 
ing 180, which corresponds to the Nyquist point fre- 
quently used in servomechanism analyses.* * When 
Yr is greater than 1/1’, for this condition, their product 
exceeds 1.0. 
the Nyquist theory in which 
transfer frequency response. It 
concluded that the 1/1), polar plot should pass over 
the origin in a counterclockwise manner with increasing 
frequency, in order to avoid an unstable airplane os- 
cillation. No attempt will be given here to explain 
the formal stability criterion, which is similar to that 


This implies an instability according to 
Y,Y, forms the loop 


may therefore be 


for the Nyquist diagram. 


(e) Calculation of Time to Damp from Combined 
Frequency Response 


Having established the condition for stability, it is 
of interest to determine, if possible, a way of predicting 
the amount of airplane oscillatory stability from these 
inverse frequency response curves. This is a compli- 
cated problem from the mathematical standpoint, es- 
pecially when all three degrees of freedom are con- 
sidered. If, however, the analysis is limited to the 
usual ranges of lateral stability coefficients of the air- 
plane, a certain correlation is possible when this system 
is reduced to an equivalent one degree of freedom 
case. From the 1/¥,4 curves generally encountered, 
the resonant frequency wz, is practically identical to the 
airplane natural frequency, which is denoted here as 
w, Also, if the damping is not too great, wz can be 
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plane with yaw damper. 


approximated by the undamped natural frequency » 
for purposes of estimating an effective static stability, 
The amount of damping should also be a function of the 
vector length R,,;,. If a damped oscillation is to occur, 
it should be expected to depend upon the airplane mass 
moment of inertia coefficient J,, the airplane damping, 
and the airplane static stability. To obtain these 
quantities, the following assumptions have been made 
to present the airplane as having an equivalent single 
degree of freedom in yaw: 

(1) The inertia will be taken as the calculated value 
of I.. 

(2) The damping will be assumed proportional to 
the ratio of the dynamic yaw response occurring at 
wp, to the effective static yaw response. 

(3) The airplane static stability giving the effective 
static yaw response can be obtained by determining an 
effective static yaw stability coefficient, Cy 5, obtained 
by considering wz as an undamped natural frequency 
We. 

The third assumption is necessary when the three 
degree of freedom case is considered. Here the airplane 
has a zero static inverse yaw response due to the neu- 
tral directional stability. This type of stability is 
neutral for all airplanes and should not be confused 
with the normal spiral stability. 

From assumptions (1) and (3), the undamped natu- 
ral frequency is 


w => QR >vV Go. a 
giving the effective static yaw stability as 


Cng = Iw,” 


The effective static yaw response then becomes 


(Wo 00) static ad —C,, c. 


NB 
since y = — £6 for the one degree of freedom case. 
From the 1/Y,4 curve the dynamic yaw response at 
WR is 


(Wo — 50) wp = 1, Rete 
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The ratio of these yaw responses will be called the dy- 
namic resonant response ratio F and is 


: (Wo — 50) wp Tw” 
sled (Yo/ — 50)static a —C,,Rmin, 

According to assumption (2), the damping should be 
afunction of F. In order to verify this, a series of sta- 
bility calculations have been made by varying several 
of the airplane stability coefficients over a wide range. 
Here, we and R»;,. were obtained from 1/Y, curves, 
and the corresponding stability equations were solved 
for the roots giving the natural frequency w, and the 
time to damp to half amplitude. In this case wz agreed 
closely with the calculated values of w,. Values of F 
were then plotted against the time to damp in periods, 
as shown in Fig. 16. In addition, a similar calculation 
was made for the yaw-damper using the simplified yaw 
damper analytical equations of Fig. 14 in conjunction 
with the airplane yaw response for a single degree of 
freedom. These results are also plotted in Fig. 16. 

These calculated values are then compared with the 
theoretical F vs. (T.,,/P) curve calculated directly for 
a single degree of freedom. In this case the airplane 


equation of motion is 


d*y . dy , , , 
i, = Cy; it + Cr = —Cx(— 4) 
when 
Cng —Cy; —C,, 
slit Al es 27 40" alia Cus 


The differential equation of motion is 


dy dy 
+ 2Far + wo” = Rswo"( — 6) 
dt? , dt v a 
where ¢ = damping factor and k; = rudder effective- 


ness. For this case it can be shown that 


On = V1 — S27 aw, wr = V1 — 27 wo 
giving 
— Dr2 - ne « / a 
a ] 2 T\/, ms 0.693 V1 — c? 
= a / at > = ‘ = 
evi-zy fF 2 i 
The airplane is critically damped when ¢ = V 1/2 = 


0.707. Any airplane that damps to half amplitude in 
a half cycle is well damped. Hence, this may be con- 
sidered a practical maximum damping. Here, ¢ = 
0.20, giving w, = 0.98, we = 0.96, or w,/wr = 1.02. 

This frequency comparison is made to show how 
close wp and w, agree for the single degree of freedom 
case. For design work with a yaw damper operating 
on the rudder, it is believed that the single degree of 
freedom curve shown, as obtained for a second order 
Stability equation, may be used to obtain (7\,,/P) for 
all lateral stability cases and still be within the neces- 
Sary engineering accuracy. 
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STABILITY CALCULATIONS FOR THE XB-47 AIRPLANE 


The method of analysis of the previous section has 
been applied to the XB-47 airplane. Here, the 1/ V4 
curves are calculated using equations of Appendix A, 
which considers three degrees of freedom of the airplane. 
These curves will then be combined with Y, curves 
obtained from laboratory tests of the yaw damper. 

It is of interest first to study the factors causing the 
1/Y, curves to change. Fig. 17 shows the change in 
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1, Y, due to different assumed degrees of airplane free- 
dom. It can be seen that the effect of roll causes the 
|/Y, curve for the three degrees of freedom case to 
move nearer the origin, thus reducing the stability. 
The lower frequency points of this curve encircle the 
origin and approach zero frequency along the +7 axis 
depending upon the sign of the é stability coefficient. 
The reason for this encirclement is the fact that at very 
low frequencies the forces due to roll are large com- 
pared to other forces present. This causes the change 
shown in the frequency response picture. 

While the low-frequency results have not been im- 
portant from the yaw damper design, they may be 
extremely important for autopilot design where the V; 
vector at zero frequency has no lead but a finite length. 
In this case the resulting 1/ Yq vector at w = 0 will be 
over to the right of the diagram, which may cause the 
1/Yqz curve to pass under the origin at very low fre- 
quencies. This may result in an airplane long-period 
instability. Fig. 17 has been drawn for IAS = 150 
m.p.h. When high air speeds are considered, the three 
degree of freedom curve can be approximated in a satis- 
factory manner for the higher frequencies by the single 
degree of freedom parabola as shown. Fig. 1Sa shows 
Y, curves where the damping in yaw, 
It may be noted 


a series of 1 
Cry» has been given different values. 
here that, if C,, j 
size, the curves of Fig. 1Sa would be shifted to the 


were increased by increasing the tail 


leit because of the C,,, change; hence, the resulting im- 
provement would be small with regard to the airplane 
damping because of the increased resonant frequency. 
The yaw damper is designed basically to increase the 
C,; without changing the frequency. Calculated 
values of F are shown on these curves to indicate the 
nature of the airplane stability, which can readily be 
noted from Fig. 16. In Fig. 18b a series of 1/ V4 curves 
are drawn for different air speeds and show how the 
higher frequency points move nearer the origin as the 
air speed is increased. These curves also show an im- 
provement in airplane stability at the higher indicated 
air speeds because J,/—C,, is reduced. The effect of 
the yaw damper servo time constant 7,, based on the 
damper frequency characteristics as defined in Fig. 
l4, is shown in Fig. 19. This figure has been drawn 
for IAS = 400 m.p.h. and shows how the airplane sta- 
bility is decreased because of increasing 7,. Next, 
the actual yaw damper frequency response will be pre- 
sented. 

The yaw damper operates the push rod, which in 
Both of these sys- 
In addi- 


turn operates the hydraulic boost. 
tems are nonlinear at the smaller amplitudes. 
tion, large airplane disturbances will cause the yaw 
damper rudder stops to be encountered at 6 = 
which introduces a large amplitude nonlinear condi- 
tion. 


5. 


Because of these nonlinear characteristics, laboratory 
Irequency response information for the entire yaw 
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damper and boost combination is necessary. These 
data are shown in Fig. 20 for two typical amplitudes. 
The yo = 
amplitude effects, which are important in studying 
The yY = +2° curve will provide 


+1/8° curve is representative of the low 


airplane hunting. 
yaw damper characteristics used to calculate the air- 
plane stability. These two polar plots represent the 
Vy feedback frequency response which will be added 
according to Fig. 15 to obtain the 1/Yj curves. In 
Fig. 21 the frequency response of the rudder boost 
alone is shown. These effects are, of course, included 
in Fig. 20. The inverse frequency response curves with 
and without the yaw damper are shown in Fig. 22 for 
two representative conditions. The 150 m.p.h. case 
represents the case where the maximum yaw damper 
stability performance is required, while the 450 m.p.h. 
condition is critical in avoiding a possible hunting con- 
These latter curves are of ititerest because of 
Of all these F 
1/8° curve 


dition. 
the large number of possible F values. 
values, the critical one occurs on the Yo = 
at we = 8.4, which shows the airplane to be stable from 
the hunting standpoint. Had a higher air speed been 
considered, this point would approach the origin causing 
an oscillatory instability. The values of F = 0.15 indi- 
cate a highly damped long period, while values of F 
near w = 4 are representative of the airplane stability 
with the damper operating. 

The original servo system developed a lagging phase 
angle near w = 3.5 and, as a result, caused the airplane 
to hunt at the higher air speeds. Fig. 23 shows the 
original servo system frequency response, together 

curves, showing that 
00 m.p.h. 


with the corresponding 1/ Y, 
hunting occurred near IAS = ¢ 


Ww pi 


EFFECT OF Bopy SIDE BENDING ON LOCATION OF YAW 
DAMPER RATE GYRO 


Recent flight tests resulted in an airplane hunting 
instability being developed at a higher frequency than 
had previously occurred (w = 14 rad. per sec.). Pilot 
observations associated this hunting with bending of the 
airplane body. The hunting condition was entirely 
corrected, however, when the rate gyro, which had been 
located near the tail, was moved forward nearer the 
airplane c.g. As a result, an analysis has been made 
to see whether this type of oscillation could be pre- 
dicted. It was assumed that the body was hinged at an 
estimated hinge location and elastically restrained such 
as to provide a body bending frequency equal to that 
found from ground vibration tests. Appendix B con- 
tains the derived equations. Fig. 24 shows the 1/V4 


and 1/Yq curves calculated for IAS = 200, 300, and 
100 m.p.h. for the rate gyro located in the tail. 
Here, the yaw damper frequency response VY,» was 


estimated by extrapolating the y = +1/8° data to 
larger frequency values. These curves indicate a hunt- 
ing near IAS 300 m.p.h., with a frequency of wp = 


17.5 rad. per sec., although this is somewhat different 


Il 
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Fic. 25. Effect of body side bending on airplane inverse fre 
quency response 1/ V4. 
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Fic. 26. Calculated rudder force and deflection during a 30 


banked steady turn. 
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from the flight record that showed that hunting hy 
occurred near IAS = 
w = 14 rad. per sec. 


350 m.p.h. and at a frequency 
The agreement is believed to be 
reasonable in light of the data used and the simple 
assumptions made. A closer analysis would requig 
a study of the actual amount of the @ motion picked 


up by the rate gyro. 


Fig. 25 shows the inverse airplane frequency 
sponse 1/ V4 to be used when the gyro is located nex 
the airplane c.g. Here, it can be seen that hunting cay 
not occur when the Y,; vectors are added to these curyes 
since the present yaw damper loses its lead near w = § 
rad. per sec. and its sensitivity is not great enough 
Had the yaw damper and control surface boost bee; 
responsive such that the Y, vector could lead until » 
18 rad. per sec., an airplane hunting could exist evey 
Here, the 
vectors, when added to the 1/ V4 vectors, could go 


for the gyro pickup near the airplane c.g. 


the opposite side and under the origin even if the 
In the light of 
these calculations, careful consideration of body bend. 
ing should be made when high responsive autopilots 
It may be con- 


sensitivity or gain ratio were small. 


are considered for elastic airplanes. 
cluded that an optimum autopilot or airplane boost 
may not be the one having the fastest response. 


YAW COMPENSATION DURING STEADY TURNING 
FLIGHT 


During a steady turn, the rate gyro deflects a con- 
stant amount due to the uniform angular velocity of the 
turn. This results in the rudder desiring to deflect 
to oppose the turn, which requires the pilot to apply an 
increased pedal force to properly position the rudder. 
The XB-47 has desirable turning characteristics when 
the damper is off, especially in smooth air where the 
The increase in the 
} 


Dutch Roll is not easily excited. 
rudder pedal forces with the damper on has been ob- 
These pedal forces have been 
While the pedal 
forces are not especially high, they give the pilots a 


jectionable to the pilots. 
calculated and are shown in Fig. 26. 


feeling of flying an airplane with a high amount of spiral 
stability, which is the actual case as indicated by the 
spiral stability boundary of Fig. 10. 

The problem of reducing or eliminating these pedal 
forces requires a means of canceling the steady or 
constant signal of the rate gyro. The basic principle of 
doing this is to provide a device that has an equal and 
opposite static signal to match the static deflection 
signal produced by the rate gyro during a steady turn. 
This device should have a slow dynamic response, char- 
acterized by a large time constant, such that it will not 
impair the desirable dynamic gyro characteristics 0¢- 
curring near the resonant frequency of the airplane. 
This time constant should be short enough, however, 
so that after a steady turn is entered the pedal forces 
will be reduced in a reasonable time. 
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LATERAL DYNAMIC 


Two ways of accomplishing this will be briefly men- 
tioned. The first method consists of including an R.C. 
(resistance capacity) electrical network in the yaw 
damper electrical system such that at very low fre- 
quencies a very small electrical signal is developed giv- 
ing a negligible rudder movement. 

A second method is to provide a slow responding 
electrical motor drive, which tends to make the electrical 
center of the gyro potentiometer follow the gyro signal. 
Here again, if the gyro is moving slowly or is positioned 
in a deflected position, as occurs in a turn, the resulting 
rudder signal will be small or zero. These methods 
are being studied, and either one should provide the 
necessary yaw compensation required during a steady 


airplane turn. 
CONCLUSIONS 


The present yaw damper has readily solved the lateral 
stability problem of the XB-47 airplane. The impor- 
tant design problem has not been one of basic operation 
of the yaw damper but rather of the proper detail de- 
sign so as to avoid phase lags occurring at small ampli- 
tudes. There are undoubtedly other ways of operating 
the control surfaces to improve the lateral stability. 
Whether these methods can be used will depend upon 
how well any new type of damper operates at small 
amplitudes. It is believed that other methods should 
be investigated with the single aim of simplifying the 
yaw damper. 

Further correlation of the present inverse frequency 
response method of analysis should be made to see 
whether the roots of the stability equation can be esti- 
mated in a similar manner for other airplane stability 
problems. 
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Appendix A 
CALCULATION OF AIRPLANE FREQUENCY RESPONSE 


IN YAW 


Lateral equations for the airplane motion are given 
These equations and notations will be 
It has been found 


in reference 1. 
used with several modifications. 
for the XB-47 airplane that the product of inertia /,, has 
a small influence on the airplane stability except at low 
speeds, where it is conservative to neglect these terms. 
For convenience, this term will be neglected. Also, 
the factor 2bu/v will be written as mp, giving the fol- 
lowing equations of motion: 


(mos — Cy,)B (—Cz)¢ (mos)P = 0 
(—C.i,)8 + (zs? — Cis) (—Cijs)\~¥ = 0 | 
7 7 , = ae” > =6(la) 
(—C,,)8 (—Crzs)o + (Is? — Cays)¥ = 
—C,,( — 6) 


From wind-tunnel tests it was found that the effect 
of the rudder motion on roll and sideslip were small; 
hence, as shown in the above equations, the rudder 
motion is considered to affect the yaw only. Ex- 
panding the left side of Eq. (1) and substituting the 
root \ for the operator s gives the airplane stability 


equation 
aod* + bod? + cod? + doAX + & = O (2a) 
where 
a = mil,I, 
bo = —(m0Ci, + Cy,Jz)I, — mol,Cy, 


Co = Cy.Ci.T. + (mCi, + Cy fr)Cr, + 
8° '@ foe) B - 
mol ,C ng — mol uC n, > (3a) 
dy =-— = Cy ,C l c ny is moC iC ng — mol i,C Ng = 
4 LC igl: = ( y gC iC ng 
& = —C1(C,.C, = G Cas) 
¥ 8 B ¥ 


The quartic stability Equation [Eq. (2)] has been most 
easily solved by the semigraphical method of Zimmer- 
man,° permitting the airplane period and time to damp 
to half amplitude to be found. 

In order to obtain the yaw response due to a rudder 
oscillation, Eq. (1) must be solved for the ratio y/—6. 
This solution gives VY, and, written in determinant 


form, becomes 





—C,, 0 
(I,s? — Cys), 0 
—CazS, Cas 
——____—_—_—_ (4a) 
—Cr, mos | 
(I,s? — Ci,5); —Cijs 


(I,s? — Cyh;s) 
¥ 


= Crs , 
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Expanding these determinants in terms of the operator 


S gives 
Vy as? + bs? + as +d, (5a) 
4=>- Va 
. s(aos4 + bos? + Cos” J+ dos + eo) 
where 
a = —C,mol, 
by = Cr,(moCi, + crt 
j Ss (6a) 
qQ = —C,,C vl lo | 
dy — Cr CiCi, / 


The frequency response is found by letting s = aj 
the resulting complex function for 
When this is done, 


and evaluating 
various values of w. 


Wo XxX, + Vij 


Y, = + “ee (Va) 
— 60 Xo + Yo] 
where 
AS = d, — byw", Y; = (0 — aya? 
Xo = — dyw? b bow’, Y> = &W — (Sa) 
Cow? + wan" 
For each frequency the vector length is 
V, = |o/—5o| = X1 cos g2/X:2 cos (9a) 
with a phase angle 
o= 1 — Pa (10a) 


where tan ¢; = ¥;/X, and tan go = Yo/X2. 

Stability coefficients for the XB-47 airplane, together 
with the resulting frequency response coefficients, are 
given in Table 1. The resulting calculations have been 
used on Figs. 17, 18, 22, and 23. 


Appendix B 


EFFECT OF Bopy SIDE BENDING ON AIRPLANE 
FREQUENCY RESPONSE 


The effect of body side bending will alter the angle of 
attack of the vertical tail. During this type of motion, 
aerodynamic loads on the airplane wing and body will 
be small compared to those developed on the vertical 
tail. For this reason the mass forces on the airplane 
body and wing will be considered, but only the aero- 
dynamic forces acting on the vertical tail will be used. 
It will also be assumed that the airplane can sideslip 
and yaw but that no roll occurs. As has been pointed 
out previously, this is a valid assumption at the higher 
frequencies. The representation of the body bending 
is made by assuming an effective body hinge pivot 
point for the aft body angular motion 0, as shown in 
Fig. 24. The after part of the body is assumed to be 
restrained from moving by an elastic restraint propor- 
tional to 6. The corresponding spring constant Ky, has 
been evaluated using the calculated weight of the entire 
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TABLE 1 


Airplane Stability Coefficients 


IAS = 150, IAS = 278, IAS = 400 

Flight Sea Level, 35,000 Ft., 20,000 Ft 
Condition Flaps Down Flaps Up Flaps Up 
V (m.p.h.) 150 500 550 
W (lbs.) 125,000 125,000 125,000 
Cr 1.5 0.443 0.214 
mo = 

2bu/V 10.4 10.1 5.36 
Iz 0.107 0.0311 0.0172 
I, 0.248 0.072 0.0346 
Cy —().1095 —(.0316 —().0289 
Ci; 0.0765 0.0107 0.0062 
Ong, —().0274 —().00198 0 
Cny —().0422 —().0112 —(). 0101 
Cig —().225 —(0).140 —(). 140 
Crs 0.16 0.12 0.12 
Crs —(0.78 —0.61 —0.61 
ao 0.276 0.0226 0. 003185 
bo 0.350 0.02788 0. OO6618 
Co 0.2726 0.0431 0.01258 
do 0.335 0.0458 0.01979 
€0 —0.00412 0.000127 0.000143 
ay 0.11138 0.0314 0. 0092 
b, 0. 1939 0.0338 0.01653 
C1 0.00854 0.00193 0.00173 
d, . 0.038375 0.0062 0.0030 


b = 116 ft. S = 1,428 sqft. 


empennage, and the measured frequency has _ been 
obtained from ground tests. 
Referring to Fig. 24, the equations of motion be- 


come: 


Side Force 


; Lf. : oe 
aoq\6 — 0 — V 0 — V y — ko) — 
W, = W _.. | oe 
rLé + VB + Vy =0 (lb 
g 4 


Yawing Moment 


, Lr. iy 
a,S,Lq\ B — 0 — 7 @ — Tad — k36}] — 


W os W a 
 L°r6 — K,2y = 0 (2b) 


Body Bending 


’ or. | a P 
a,SqLr \ B — 6 — y @ — yp — k36}) — 


W . W W a 
* Lr) — — Lf + — LrVB + 
g g g 
W, - . 
LrVi —K,@=0 (3b) 
g 
where 
a, = slope of vertical tail lift curve (C, per rad.) 
S, = vertical tail area, sq.ft. 
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L = tail length 

WW = gross weight of entire airplane, Ibs. 

W, = weight of empennage, Ibs. 

K, = radius of gyration of complete airplane 
K, = body spring constant, ft.lbs. per rad. 


These equations will be written in a nondimensional 
form for convenience of solution. Basing all lengths 


in terms of half wing spans, the following notation is 


defined : 
m = 2(W/g)/pSb 
a = 2 W./g)/pSb 
k, = 2K,/b 
Crug = 40Si/S)(L/6) 
Cy, = —a,(S,/S) 
m = Cy, ‘Qu 
E = K,/qSbr 
l = (2L/b) 
k; = —Cn, Cn, 
De = (2V/b)6, D0 = (2V/b)*6 
De = (b/2V)6 


Substituting the above notation in the equations of 


motion and rearranging, 


(D8 — mB] + [(1 + ml)Dy) + 


[- ** rD* + mlrD@ + m6] = mk;(—6) (4b) 
m 
[Cag] + (wh2DY + CyglD¥] + 
[ul2rD° + CydrD0 + Cr] = Cugks(—8) (5b) 
[-udDB — C,,8) + [ud?D yy + 
(Cr, — widlDy) + [ud?rD°@ + C, JDO + 
(E+ Cn.)4| = Cr Ral — 8) (6b) 


Based on the above equations, the stability equa- 


tion becomes: 


dy’d4 + bo’A? + co’? + dh’A + & = O (7b) 
where 
%' = wrk? (uw — wr) — ple 
bh’ = C,, Jr [k2?(u — pw) — pdl*] 
a = Cx ie r (1 1 #, 
. E F Mr r + (8b) 


miu (pk? — wl?r) + pkR2E 
d)’ = Cn LE — mpkek 
@ = Cr.E 


p 


Eqs. (7) and (8) permit the periods and time to damp 
to be calculated using the following equations: 
h=G+ aj, w= (2V/b)a (oF 

y , i veu = ant JD) 
P= (2r/a@)(b/2V), (—0.693 /a)(b/2V)§ 

Equations giving the inverse frequency response can 
be calculated by the procedure given in Appendix A. 
Here, 
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l —5o [e0’ — co’@? + ao’w*|] + [do’a — do’a*)j 
V4 Wo cy" 
(10b) 
where ¢;’ = RsCr,E. 
The body bending inverse response is 

— do [eo’ = Co’ @? a dy’ @*} + [dy’@ —_ bya | 7 i 

” = (11b) 
A —Ay'w? 
where 


dy! = pk-ks(mpd + Cn,) ss HiCn Rel? 


It is desired to calculate the effect of locating the yaw 
damper gyro in the tail where the gyro pickup is 
(0+ ); hence, the airplane inverse frequency response 
to be used for the gyro 1n the tail is 


oS ae <te 
V4 tail 7 (90 T Yo) 


[@o’ — co’@? + dy’) + [do’@ — bo’w*)]j 


[c,’ we Ay’ wo? | 


(12b) 


Eq. (10) is used to calculate the 1/4 curves of Fig. 
25 and Eq. (12) gives the 1/ Y, curve for the rate gyro 
located in the tail portion of the airplane as used in 
Fig. 24. 

Typical coefficients for the XB-47 airplane for IAS = 
400 m.p.h. at 20,000 ft. are given below: 


W = 125,000 Ibs. W, 8,477 

b = 216%. r = 0.333 

S = 1,428 sq.ft. wo = 138 cycles per min. 
= 14.4 rad. per 
sec. 

p = 0.001265 S, = 227 sq.ft. 

KK. = 246K. (b/2V) = 0.072 

L = 46ft. l = 0.792 

Then 
m 37.0 m = —0.00762 
% = 2.51 Ky = wo(rL)?(W,/g) = 
12.8 *& 10° 

Cy, = —0.564 E = 0.568 

Crg = 0.228 k,? = 0.181 

kg. = 0:625 


giving 


As’ = 0.6044 
Cy’ = 0.0734 


Qo’ = 2.448 


bo’ = 0.2747 


rd = 1.937 
dy’ = 0.1293 
é&’ = 0.1267 


(Continued on next page 
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Appendix C 


SERVOMECHANISM DERIVATION OF 
INVERSE TRANSFER METHOD 


In this appendix the method of obtaining 1/ Vz from 
the servomechanism standpoint will be analyzed. Fig. 
6 is a block diagram of the airplane rudder system with 
hydraulic boost and yaw damper. 

The only physical input to this system is through the 
pilot’s pedal motion (—ép). The output motion is the 
yaw y. The feedback signal 6, opposes the pilot pedal 
motion; hence, the rudder boost output —6 for the 
closed loop is 


—6 = Yz(—6- — 6p) (Ie) 
The various transfer functions may be expressed as 
follows: 
Y = Y,4(—64) airplane transfer function 
—6 = Y,(—6p) rudder boost (no feedback included) 
by = Yr(y) feedback transfer function 


When these relations are substituted in Eq. (lc) for 


the closed loop circuit, the output y due to the input 
— 6 may be written as 

output y V,Yp ; 

input —dp 1+ Y,YsY, " 

It is desired, however, to consider the rudder motion 

as the input and to modify this by the response of the 
feedback to form 1/Y 3. As the rudder can only be 
operated through the boost, a signal must be applied 
to the rudder pedal to produce the reference input 
—6 considering the feedback 6, absent. This input 
pedal motion is —ép = —6/Y,; then the yaw response 
with the feedback operating becomes 


¥/-6 = Vg = Vu/( + Ya¥o¥r) (3¢ 

or the inverse airplane frequency response with yay 
damper operating is 

1/VYz = (1/¥a) + YVrY;z (4c 


In the analysis considered in this report, the hydraulic 
boost transfer function ), has been included with the 
feedback response, and the product Y,J>y has been 
used as Vy. 
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Steady Flow of Nonvi 


scous Elastic Fluids in 


Axially Svmmetric Channels 


M. H. VAVRA* 
United Slates Naval Postgraduate School 


SUMMARY 


Subject matter is the axially symmetric steady flow of a non- 
yiscous elastic fluid, with constant totai energy, between two 
wbitrary coaxial surfaces of revolution. By means of orthogo- 


gal, curvilinear coordinates, a transformation of the Eulerian 
equation of motion is obtained which permits the determination 
of the meridional velocities for arbitrary distributions of the 
tangential velocity components. Based on this relation, an 
iteration method is described which serves to calculate arbitrary 
three-dimensional flow fields 

An application of the theory to an irrotational flow in an inlet 
channel for a centrifugal compressor wheel shows the importance 
of determining the flow ahead of rotating wheeis in order to 
evaluate the correct inlet conditions for the impeller 

Rotational motion gives flow patterns that are entirely differ- 
ent from those of irrotational motions. Depending on the vor- 
ticity and the shape of the channel contours, there exist definite 
relations between the meridional and the tangential velocity 
components to make possible a stable motion that follows the 
hannel surfaces without separation. 

The method is particularly suited to determine pressure gra- 
lients in curved diffusors and inlet ducts of compressors and tur- 
bines in order to estimate the possibility of flow separation 


SYMBOLS 


E = total energy per unit mass 
mass flow per unit time 


S dp p 

V = magnitude of velocity 

V; = magnitude of tangential velocity 

V2 = magnitude of meridional velocity 

a = velocity of sound at stagnation 

hy, he, h = variable parameters of the system of coordinates 

n = distance along normals to the meridian lines 

p = static pressure 

pb = stagnation pressure 

r = distance from axis of symmetry 

R = radius of curvature of meridian 

li, U2, 43 = orthogonal coordinates 

J = velocity vector 

é = vector of external force per unit mass 

i, Io, 1; = unit vectors in directions 1, M2, “3, respectively 

j = unit vector in direction of the axis of symmetry 

f = radius vector 

. = potential of fe 

2, = scalar, component of V X V in direction 1 

thy = scalar, component of V X V in direction “2 

Y = Cp/Cy = ratio of specific heats at constant pressure and 
constant volume 

0 = angle of surface u, = constant with fixed plane, 
containing the axis of symmetry 

p = mass density 

po = mass density at stagnation 


Received February 1, 1949. Revised and received July 15, 


1949, 


* Professor, Department of Aeronautics 


Subscripts 
1 = tangential direction 
2 = meridional direction 
3 = direction of normal 


INTRODUCTION 


A NONVISCOUS ELASTIC FLUID is made to undergo a 
steady motion in an axially symmetric channel 
which is formed by the rotation around an axis of the 
curves Cy and C; of Fig. 1. The channel is fixed in 
space and does not contain blades to guide the flow. 
The total energy of each fluid particle is the same and 
no heat exchange shall take place between the fluid 
the channel contours. The meridian line m of 
rotated around the axis, produces a surface of 


and 
Fig. 1, 
revolution to which the flow velocities shall be tangent. 
The curves, perpendicular to the meridians, are called 


normals. 
DERIVATION OF FLOW RELATIONS 


Equation of Motion 
The Eulerian equation of motion for steady flow is 
from reference 1, page 62, 
V-VV = fe — l/p)Vp (1) 
or 


V-VV = -V(P+Yy) (2) 


since the only forces acting upon the mass particles are 
those due to gravity. They have a potential; hence, 


) 


fe = -Vy (3) 


With P = /dp/p, there is 


(1/p)Vp = VP (4) 


If V denotes the magnitude of 7, Eq. (2) may be ex- 


pressed by 
VxvxV=v[(V2/2)+ P+) =VE (5) 


E represents the total energy per unit mass, consisting 
of the kinetic energy, the enthalpy, and the geopoten- 
tial energy. Since the total energy is constant at all 
stations in the flow, there is VE = 0, which implies 
that the flow must either be irrotational or that the 
curl 7 must be parallel to the velocity vector V at all 
stations in the flow.” 
From Eq. (5), for VE = 0, and Eq. (2) 


V-VV = V(V?/2) (6) 
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i, 72, 73 are the unit vectors tangent to the coordinate | fron 
U1, U2, U3. They forma right-handed system, as show 
in Figs. 2 and 3. From these figures, there are, sing 
u, = 06, ince 
syiml 
hy = PP, dit, = rdé dé = r °. 
to tl 
hy = PP2/duz = ds/dus 10 Fig. 
h. = PP; dus = dn dus 
Equation of Motion in Orthogonal Coordinates Fron 
There are 
VP =4Vi+ nl? D —_ 
Simi 
V2 = V+ Ve? (13 
All partial derivatives of scalar values, with respect 
to m1, are equal to zero because of the axial symmetry oj 
the flow. With the operator V from reference 1, pag 
15, ; 
Nt 
so V1? 07 ViVe0%m = WiV20% and 
hy On, hy Ou, hy Ou ds\ ) 
V2? Ot. _ hs _ 1 O(V2?/2) ) 
P : TT UY ‘ 12 : - 
Fic. 1. Meridional contours of axially symmetric channel hy Ons he Our hy Ou, iV, 
with arbitrary station P on meridian m. 
and 
~ ‘ ; , ; , ; : 2 ; 72/9 ms 2/9 
Eq. (6) describes the motion of a fluid particle in a flow v(_ ) _ 2 o(Vi*/ 2) 4 is O( V2*/2) 
which has constant total energy. Moreover, the flow 2 hy Ota hy Ole 
must be steady, and, similar to the assumptions made is O(V1?/2) 73 O( V2?/2) Z 
for the derivation of the Eulerian equation, the vari hs us hs us The 
ables and their derivatives must be single-valued, finite hai 
and continuous functions in the flow field under con- The changes of the unit vectors in direction of the 11) 
. . . ~ . (g\ 
sideration. Therefore, Eq. (6) cannot be applied to a ordinates can easily be determined. For instance, won 
region in which there exists a shock front. It is, how- the increase of 7 from m to m + di equals I,( ) 
. e . ° . . 3 
ever, applicable for elastic fluids independent of whether diy = (O%/0m)duy lf to th 
the flow is rotational or irrotational. f 
y 
Orthogonal Coordinates ee 
. ; a" 7. Ot 
An arbitrary station in the flow region—say, the _ i+ 3 dus 
point P of Fig. 1—shall be located by the three mu- C Ous rt 
tually perpendicular, coordinate surfaces uw; = con- betw 
stant, #2 = constant, and uw; = constant, as shown in lines 
Fig. 2. The coordinate surface uz = constant is the curv 
surface of revolution, obtained by rotating the meridian the 1 
m around the axis of symmetry. The surface uw. = Rat 
constant is a surface of revolution with a generatrix Aen 
perpendicular to the meridian lines. Fig. 3 shows that _— 
the coordinate surface “4; = constant is a plane contain- in Fi 
ing the axis with w; = @ = constant. Ec 
Let ¢ be the radius vector of the point P from a fixed tions 
origin 0 on the axis. A neighboring point to P shall be 
located at the terminus of § + dé. Thus, 
— o¢ o¢ 
dé = —- du, + — du. + = dug (7) f 
Ou; Ou> Ou; 
; Ei 
. : / + = ° 
or with the variable parameters /1, ho, hs l Uz Cuz Cons t. dion 
dg = ihydu, + igheduz + ishsdus (3) Fic. 2. Orthogonal coordinates in meridian plane. tivel 
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From Fig. 3 


di, = jd0 x 1; (17) 
since the unit vector 7, is rotated around the axis of 
symmetry by an angle d@. The unit vector j is parallel 
to the axis of symmetry, with the direction shown in 


Fig. 2. There is 7 = 72 sin \ + 73 cos \, and 
di, = d0(sin X\ 72 K % + cos XA 23 XK 1h) (18) 
From Eqs. (16) and (18), since du, = dé, 
07,;/0u, = —i3 sin X + 7 cos A (19) 
Similarly, 
Oi2/OU2 = 73 (dX/du») (20) 
07/0. = 0 (21) 
Oi2/OU; = —1, COS A (22) 


With the variable parameters of Eqs. (9) to (11) 
and writing [O( )/Ow2] du, = do( ), [0( )/Ous] dus = 
d;( ), Eq. (6) becomes 


: V1 do( V1) pee V1 
iV) — cos \ + + DV; cos \ — 
r ds r 
ds a f | Vi? — y,20 
re + 13 ; sin T 2 a 


ds( Vi? 2) d( V2" 2) 
— =@Q (23) 
dn dn 


The differential d.(_) represents the increase of a scalar 
from “2 to uw + dus or along ds in Fig. 2. Hence, 
d,( )/ds equals the rate of change of the scalar with 
respect to the meridian Similarly, 
d;( )/dn is the rate of change of a scalar with respect 
to the normals to the meridian lines or + 173. 


+ uo. 


lines or 


From Fig. 2, 


dd\/ds = —1/R (24) 

The minus sign is introduced to make a distinction 
between convex and concave shapes of the meridian 
lines. A radius vector, with its origin at the center of 
curvature C and its terminus at the point of reference at 
the meridian, is introduced. The radius of curvature 
Rat a particular point on a meridian line is defined as 
a positive value if the aforementioned radius vector 
points in direction of +7; and vice versa. Hence, R 
in Fig. 2is a negative radius of curvature. 


Eq. (23), with Eq. (24), establishes the following rela- 


tions: 
(V,/r) cos X = do(V;)/ds (25) 
ait. Sl d3( V1") + 9 V2" ‘ d3( V2") 0 (26 
4 sin t Z = (20) 
r dn R dn 


Eqs. (25) and (26) are the equations of motion in 
directions of the meridians and the normals, respec- 
tively, which can easily be integrated. 
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Ui+dui=const. 


i1,2bdu, i 





plane perpendicular to 


coordinates in 
the axis of symmetry seen from the origin 0 toward A of Fig. 2. 


Fic. 3. Orthogonal 


Integration of the Equations of Motion 


From Fig. 2, 


ds cos \ = —d,(r) (27) 
Eq. (25), with Eq. (27), 
rds( Vi) + Vyds(r) = ds( Vir) = 0 (28) 


Eq. (28) states that the angular momentum per unit 
mass remains constant along uw or on a surface of revolu- 
tion produced by the rotation of a particular meridian 
m around the axis of symmetry. This relation is ob- 
vious, since no external moment is exerted upon the 


fluid particles. From Fig. 2, 


sin \ = d3(r)/dn (29) 
Hence, Eq. (26) may be expressed as 
ds( Vo") Z V; ds( | 1") ? 
| 2~ : Q (30) 
dn R r dn 


This linear differential equation of the first order for 
V2, with variable coefficients and the variable n, has 


the general solution‘ 


aii ; V; ds( Vir) . J 
R)dn} CH / 9 e + JS (2/R) dr dn 
. r dn 


(31) 


V. 9 a 


9 = €é 


Therefore, it is possible to determine the distribution of 
the meridional velocity component | along the nor- 
mals if r, R, and V; are known. C is a constant of 
integration. 

Eq. (31) can be modified by introducing the vorticity 
vector V X V of the flow. From reference 1, pages 14 
and 16, there is for the variable parameters of Eqs. (9) 
to (11), 


with 
Vi 1 ds| Vir) 


Vor dn 
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normal n* 

















and 


QQ = 


(1/r) [d3( Vir) /dn] (34) 


The curl of V has no component in direction of the 
normals, and, since 2)/Q, = V;/ Ve, it can be seen that 
the vector V X VP is always parallel with the velocity 
vector /. 

Irrotational motion exists if either ; = O or 
d3(Vir)/dn = 0. In the first case, the flow has no tangen- 
tial velocity components. The other case represents a 
flow for which the product Vir is constant along the 
normals, or, since |\r is constant anyway along the 
meridians, a free-vortex flow with V;r = constant at all 
stations in the flow. The meridional component Q of 


V X V may be substituted into Eq. (31). 


Distribution of the Meridian Velocity V, Along the Normals 


For particular channel contours, there shall be known 
the meridian lines (or surfaces u3 = constant) and the 
normals to the meridians (or surfaces “2 = constant). 
As shown in Fig. 4, the vorticity shall be given at one 
The vector V X VP at station 
P on an arbitrary normal n can be determined from 
V X V* at the point P* on the normal n*, and on the 
same surface uz = constant as the point P. 


particular normal n*. 


Applying 
Gauss’ divergence theorem (reference 1, page 20) for 
V X PV toa region bounded by the surfaces us = con- 
stant, uw; + du; = constant, and the surfaces uw. = 
constant at the normals m and n*, gives 


since V-V X V = 0. The vector da in the surface 
integral is a normal surface vector of an element of 
area and is considered positive if it points away from 
the bounded region. With Eq. (32) and 
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Eq. (35) gives, with the notations of Fig. 4, 
da-V X V+ da*-V X T* =0 
or, in scalar form, 
Q. = (r*/r) (dn*/dn) Q* (37 


Hence, it is sufficient to know the vorticity along a par. 
ticular normal n* and the shape of the meridian ling 
to obtain the vorticity at arbitrary stations in the floy 
field. 

From Eq. (28), 


Vy" (r*/7) (38 


V;* is the tangential velocity component on the normal 
n* at the point on the same surface “3 = constant as th. 
station where the velocity V; exists. 
The velocity V2 shall be determined at the station / 
Meridian lines and normals are known 
V,* and, consequently, 22* 
Let Q be a moving point 


of Fig. 5. 
Moreover, are given along 
the particular normal n*. 
on the normal n between the surfaces “3; = constant 
and wz; = constant. To P and Q correspond the 
points P* and Q* on the normal n*. 
of Eq. (81) shall be carried out in direction +73. Thus 


the constant of integration of Eq. (31) is 
C = ( Vo; iz 


namely, the square of the meridian velocity on the nor 

mal at the surface “3; = constant. 

of Fig. 5, Eq. (31) becomes, with Eqs. (34) and (38), 
V2? = A*(V2,? — B?) (39 


wherein 


At=e “"=' (40 


The integrations 


Using the symbol 
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an* *!\9 n’ 
anes : (7 ‘= oSu —,@n R 
p= 2f V."'0,* —— dn* (41) 
n* =0 ae i 


,2*’, r*’, n* are values corresponding to the point 
0* on the normal n*, whereas 7’ and n’ correspond to 
station Q on the normal nm. The first integral in Eq. 
41) can always be evaluated along the normal n* from 
= constant to uw; = constant, independent of the 
However, the 


Msi 
normal where the point P is located. 


integrals constituting the powers of the basis e must be 
slved for the particular normal » on which P is situ- 
ated. 

Eq. (39) holds for flow fields that satisfy the same 
conditions as those necessary for the use of Eq. (6). It 
is interesting to note that neither the density of the 
fuid nor the potential of the gravity forces appears in 
fq. (39). Fora given pattern of the meridian lines and 
the distribution of the ‘tangential velocity components 
along one particular normal, Eq. (39) makes possible 
the evaluation of the meridional velocity at all stations 
in the flow, provided that the velocities V2; along the 
inner contour are known. 

For irrotational flows—i.e., for free-vortex flows or 
for V; = O—there is 2, = 0, and, from Eq. (39), 

™n /p 
Ve/Van = @ Se (42) 
This relation has been derived by G. Fluegel in 1914.° 
Another particular solution is obtained from Eq. (45) 


if, besides Q2 = O, the centers of curvature of all merid- 
ian lines passing across a normal coincide. In this 
case, the normal is a straight line, and 
dn = dR 
or [see reference 6, Eq. (329), page 291] 
VoR = Vo,R, = constant (43) 


R, is the radius of curvature of the inner contour of 
the channel. For the aforementioned conditions, not 
oly the angular momentum of V; with respect to the 
axis of symmetry is constant but also the angular 
momentum of V2 with respect to the center of curva- 
ture. 

The velocities V2; can be determined with the equa- 
The mass flow across a normal is 


n 
on f pr Vodn 
0 


y M/29r 
- Soper AV 1 —(B/V24)2dn 


tion of continuity. 


M = (44) 


and with Eq. (39) 


(45) 


for an incompressible fluid there is p = po = constant. 
For an isentropic process of a perfect gas, neglecting 
seopotential energy differences,° 
9 r \ FI/(y-1) 
p y¥—-1/V;? V,? ” 
es 2 a ae - + — (46) 


9 2 2 
Po a ao ao 
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Eqs. (45) and (46) show that the velocity V2; can only 
be determined directly if B = 0 and p = pp—i.e., for an 
irrotational flow of an incompressible fluid. For a 
general case, a value of V2, must be chosen to deter- 
mine V2 and p from Eqs. (39) and (46), respectively. 
Thereafter, the chosen V2; has to be checked by means 
of Eq. (45). 


CALCULATING PROCEDURE 


An iteration method must be used to determine the 
flow field in a given channel, since the meridian lines of 
the flow must be known before the above relations can 
be used. First, an approximate picture of the meridian 
lines must be established, and it is useful to apply meth- 
ods that are indicated on pages 136 and 291 of reference 
6. However, in most cases it is possible to anticipate 
the ratios V2/ V2; along the normals from the curvatures 
of the boundaries and the component Q) of the curl 
vector by properly interpreting Eq. (39). The merid- 
ian lines of this first approximation are preferably 
drawn in such a manner that equal mass flow exists 
between two neighboring surfaces “u; = constant. 
Further, the radii of curvature of these meridians have 
to be determined along different normals. If graphical 
methods are used, this procedure has to be carried out 
with great accuracy to obtain consistent results. 

From the flow conditions at the normal n*, it is then 
possible to determine the values of A and B, with Eqs. 
(40) and (41), for the first approximation of the flow 
lines. Subsequently, the values V2, of the meridian 
velocities at the inner contour must be chosen for the 
different normals 
Eq. (45) is used to determine whether the value of 
2;, necessary for the given mass flow through the 
channel, equals the initially chosen value of V2;. The 
assumed magnitude of V2; must be changed until this 
agreement is obtained. With these data it becomes 
possible to determine the points of intersection of the 
meridian lines with the assumed normal for equal mass 


J re 


flow between surfaces vu; = constant. 

This procedure, carried out for different normals, 
gives a new set of flow lines which will probably differ 
from the first approximation. The new flow lines will 
then be used as the second approximation with which 
the afore-described procedure has to be repeated. In 
this manner successive approximations have to be car 
ried out until the assumed flow net coincides with the 


calculated one. 


NUMERICAL EXAMPLE OF AN IRROTATIONAL FLowTt 


Fig. 6 shows an inlet channel ahead of a centrifugal 
compressor wheel with a weight flow of 35 Ibs. of air 
per sec. Stagnation pressure and temperature of the 


+ The calculations of this example and the mathematical 
analysis of the channel contours were carried out by Lt. Gerald 
M. Monroe, U.S.N., now Graduate Student, California Institute 
of Technology. 
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2 1/2 in. s = 60 a, 0 
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j-—<— From the two points at stations A and J, it is pogsibj ~ 
y4 Lut I ~ —— to determine the fixed origins 0 and the angles 6,, where 
TTN EF a i the curves have zero curvature. For the inner conto geth 
U30 _ 4 | t or the surface 73; = constant, there is rota 
yy el ust 
U3 s = (11.84) 6-167) with 6, = 21.4 nl 
u3” —3 . 
; > re cale 
Us? — For the outer contour 
sec. 
~ = (997 (0.0402) wos = : 
“ s (9.27) 0 with 6; 11.4 diffe 
geen . in. ‘a ° “s . 7 ¢ 
Us3i , From these equations, the radii of curvature of the by : 
contours can be calculated. illus 
— ; upol 
Vir is constant at all stations and equals 46.6 sq.ft. per : 
. . ° an¢ 
sec. Preliminary calculations showed that p/p» fron th 
: E Y F ec 
Eq. (46) varies less than 3 per cent between J-J and 4. The 
A. Therefore, constant density was assumed along the 
particular normals, and the respective ratios p/» 3 of 
were determined from the mean values at each nor ss . 
. 
mal. 
aXxla 
“7 — For the first approximation of the meridian lines, it cons 
Fic. 6. Flow lines in channel with irrotational motion. was assumed that the curvatures of the meridians at angl 
station A-A change linearly from 1/R = (0.0932) in onal 
, aT ; wee atr = 4in. to1/R = (0.1097) atr = Yin. At station Ir 
air are 14.7 lbs. per sq.in. absolute < F., respec- a: 
ur are 14.4 bs. per sq.in bs lute ind 70°! * Fespec- — J-J the meridian lines are parallel and equally spaced aes 
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components along station J-/. The curv: ‘s of the PF ‘ «of 

. I on along station I he wy itures of the meridian lines that satisfy the relation s = bO* be J nig 

channel contours are zero at station J-J : rradually . e 
anne! contours are zero at station ind gradually tween the two known points at A-A and J/-J, respec 
















increase to a maximum at station A-A. Such con- tively, with zero curvature at J-J. By drawing th i 

engl may be expressed by equations of the form (see jormals with the calculated slopes and using Eq. (42 helo 
Fig. 6) with the determined radii of curvature the ratios 
V2/ V2; along the different normals are given. Conse 
1.0 : ‘ quently, the points of intersection of the second approxi 

mation of the meridian lines with the initial normals can T 

be established, again for equal mass flow between sur- st s 

faces v3 = constant. These points are marked by the that 

0.95 = | | | L_| small circles in Fig. 6. It can be seen that they prac arbi: 

f tically coincide with the corresponding points of the with 
P/p, 4 first approximation. Thus, no repetition of the calew- 

| = — DM oe lating procedure seems warranted. 

0.90 -20— , eae tt In Fig. 7 are plotted the pressure ratios // po along The 

V2/v2* es the inner and outer contours of the channel. po is the stat 

| } stagnation pressure of the fluid. The pressure gradient of a 

[ at the outer contour reaches values that are about tour 
=n / TI times the maximum pressure gradient at the inner con “a 
| | Pi | tour of the channel. Fig. 7 also shows the velocity aia 
| ratios V2/V2* along the contours. The constant : 
1 nee es Ce 2, meridian velocity I,* at station J-/J equals 218 ft. per qT 
T m . Sup] 
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J H G F E D C B A V. t st ti { 1 f 60.5 ft t the lines 
i » at station A-A changes from 260.0 It. per sec. at : 
distance along contour/|—- te ang i eae 
radius r = 4 in. to 426 ft. per sec. atry = Qin. Fig. 
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Fic. 7. Velocities and pressures along the contours of the : Mor 


channel of Fig. 6 for irrotational motion. k-E, and J-J of Fig. 6, together with the flow angles 
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a, measured between the actual velocities and the cir- 
cumferential direction. 

In Fig. 9 are shown the angles a at station A-A, to- 
gether with the flow angles 8, relative to an impeller 
rotating at 6,000 r.p.m. The angles a and £6 are 
ysually determined for constant meridional velocity at 
the impeller inlet.’ In Fig. 9, ag, and By, are angles 
calculated with a constant velocity V2 of 345 ft. per 
sec. It can be seen that the actual relative flow angles 
differ from B,y. by about 8° at the inner radius and 
by about 6° at the outer radius. Fig. 6 of reference 7 
illustrates the considerable influence of the angles 6 
upon the efficiency of centrifugal compressor wheels, 
and it seems necessary that proper attention is paid to 
the actual velocity distributions at the impeller inlet. 
The ‘angles of attack,”’ proposed in reference 7, are of 
the order of magnitude of the difference between 8,,, and 
sof Fig. 9, although the present inlet channel has ex- 
tremely small curvatures. Usual designs, having small 
axial extensions and bigger curvatures, will produce 
considerably larger differences between the actual 
angles 8 and the ones calculated with constant meridi- 
onal velocity. 

Inlet ducts, similar to the one of Fig. 6, produce 
meridian velocities at station A-A which increase with 
increasing radius. Especially with channels having 
small axial extensions, the actual Mach Number of the 
relative flow at the outer diameter of the impeller inlet 
will be considerably higher than the one calculated with 
an assumed constant meridian velocity. This fact 
might explain the phenomenon that ‘“‘choking”’ effects 
are sometimes experienced with compressor wheels 
having apparent local Mach Numbers that are well 
below unity. 


NUMERICAL EXAMPLE OF A ROTATIONAL FLow 


The channel of Fig. 6 shall be used to produce a flow 
at station A-A which has tangential components |; 
that are proportional to the radius r. With ¢ as an 
arbitrary constant, there is 1; = cr, and from Eq. (34) 
with dn = dr 


The meridional component of the vorticity vector at 
station A-A is therefore constant and equal to the curl 
of a flow which rotates around the axis in the same 
manner as a solid body. Station A-A will be used as 
the normal * along which the integrations of Eq. (41) 
are carried out. 

The meridian lines of the first approximation are 
supposed to have horizontal tangents at station A-A 
and vertical tangents at J-/ of Fig. 6. The meridian 
lines at J-J will have certain curvatures that, however, 
are small since the contours are straight lines, and the 
integrals of Eq. (40) will very nearly equal unity. 
Moreover, the assumption is made that the curvatures 
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Fic. 8. Velocity distribution and flow angles, at the normals 


J-J, E-E, and A-A of the channel of Fig. 6, for irrotational mo- 
tion. 


along A-A change in the same manner as those of the 
calculated meridian lines of the foregoing example. 
With a weight flow of 14 Ibs. of air per sec. of the 
channel and the same air conditions as in the previous 
example, the density changes may be ignored. A and 
B along station A-A can be determined with Eqs. (40) 
and (41). With a trial-and-error procedure it will 
then be possible to evaluate the velocity V2; from 
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Fic. 9. Actual flow angles a and 8 for calculated velocity dis- 
tribution at station A-A of Fig. 6, for irrotational motion and an 
impeller at 6,000 r.p.m., and angles ayy, and Bay, for constant 
meridian velocity. 
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Eq. (45) and, hence, also V: with Eq. (39). For a 
chosen value of 1; = (127.32)r (Vi in ft. per sec., 7 in 
ft.) at station A-A, there was V2; = 125 ft. per sec. 
Thus, the points of intersection of different meridians 
with station A-A can be established. Eq. (41) can 
now be solved for station J-/, with r* and dn* from 
station A-A. Trial-and-error methods, with chosen 
values of V2;, introduced in Eq. (45), give the points of 
intersection at this station along which the radius 7 is 
constant and R equals infinity by assumption. 

The first approximation of the flow net shall be again 
a set of meridian lines between the two known points at 
stations A-A and J-J satisfying the relation s = 
bé-*. With the radii 7 and R of these flow lines, it is 
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possible to establish a second approximation in the Ve/Vei | | 
manner outlined before. Also in this example, a good Ui?” Uso Usi >) «Uso Uni Us’ Us” Us" Uy 
coincidence between these two flow nets was reached, 0 | = 
making further approximations unnecessary. —j-i— ee | le A-A sal 


Fig. 10 shows the velocity profiles and the flow angles : ey ee ee 

‘ i a ss Fic. 10. Velocity distribution and flow angles at the normals 
a at stations A-A and J-J of Fig. 6, as well as along a = J-J, E-E, and A-A of the channel of Fig. 6, for rotational motio 
normal almost at the same position as E-E in Fig. 6. With constant vorticity at station 4-4. 
The meridional velocities V2 are practically constant 
along the three stations, in spite of the varying radii of the velocities will not be able to produce the forces 

‘“ = awe > » 6 > “4° . . 

curvature. Contrary to the first example, the angles yecessary for the equilibrium of the mass particles 
a are decreasing from the inner to the outer contour The flow would separate from the outer contour of the 


of the channel. channel, thereby increasing the meridional velocity 
A peculiar condition can be seen from Eq. (39), since components until equilibrium becomes again possible. 
4 D4 y = © $ D4 « > © a eC . . . *,* 

V2 equals zero if V2; = B and, if B is larger than V2, For free-vortex flows, no such limiting conditions 


equilibrium of the forces acting upon a flow particle is components can be increased arbitrarily. In actual 
possible. Since % and, consequently, the value of B jrrotational flows, however, where the relation Vr = 
from Eq. (41) depends on the distribution of the tangen- constant does not hold within the boundary layers, 
tial velocity component along the normals, it is not pos- there seems to exist the possibility of a minimum flow 
sible to choose the magnitude of V; arbitrarily for a rate, below which the flow might become unstable. 
given massflow through the channel. Toeach mass flow Sycha phenomenon could be the reason for flow separa- 
and a given vorticity at a particular station, there cor- tjons that are different from those due to pressurt 
responds a maximum value of V;, for which the meridi- gradients at the boundaries. 

onal velocity component becomes zero at one or more 

stations in the flow. Conversely, for a given V; dis- GENERAL REMARKS 

tribution, there exists a minimum mass flow that pro- 
duces the same effect. In the present case, the maxi- 
mum values of B exist at the outer contour. There are: 


the values of V2 are imaginary. This indicates that no exist, since B always equals zero, and the tangential 


Considerable calculating work and numerous graphi- 
cal integrations are required to carry out the proposed 
method. Especially for channels having large changes 


at A-A: B= 94.2 ft. persec. V2, = 125 ft. per. of the curvatures of the contours, hence, where it is not 
a peree . _ _ 5€€+ possible to select a reasonably close first approximation, 

at E-E: B= 76.6 ft. persec. V2; = 108.7 ft. a number of successive trials will be necessary. How- 
P : per sec. “ever, it was found that the results converge rapidly for 

atJ-J: B= 42.7 ft. persec. Vo = 90.4 ft. flows at low Mach Numbers. Supersonic flow im 


per sec. axially symmetric channels could be treated with the 
method set forth, provided no shock fronts exist, which 
will roughly be the case if the meridian velocities are 
subsonic for supersonic total velocities. However, itis 
possible that the method converges but slowly, or not 
at all, if the first approximation of the meridian lines 


Hence, with the same tangential velocity components, 
it is only possible to reduce the flow rate from its as- 
sumed value of 14 Ibs. per sec. to about (14) (94.2)/125 
= 10.5lbs. per sec., for a reduction to zero of the merid- 
ian velocity at the outer contour at station A-A. At 
E-E and J-J, the smallest possible flow rates would be differs too much from the actual flow net. 

Despite these limitations, the method will give inter- 


9.95 and 6.65 lbs. per sec., respectively. Therefore, 
esting results for many applications, such as diffusors 


conditions at A-A limit the minimum flow rate. For 
still smaller mass flows the flow will be unstable, since (Continued on page 172) 
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On a Family of Oblique Shock Curves 


EDWARD R. C. MILES* 


Institute for Cooperative Resear 


SUMMARY 


The curves .V sin o = constant, where / is the free-stream 
Mach Number and o is the shock angle, are level curves for the 
principal pressure, density, and temperature ratios across an 
oblique shock. If these curves are adjoined to the family of 
shock polars, dynamic relations, as well as the kinematic rela- 
tions associated with the shock polar, are exhibited graphically. 
It has been customary to obtain these curves by interpolation. 


In this note, a simple analytic representation is derived 


INTRODUCTION 


I" IS KNOWN that the principal ratios of state vari- 
ables behind and in front of an oblique shock are 
functions of 7 sin o, where .V/ is the free-stream Mach 


Number and o is the shock angle.’ * Thus, with 
7 = 1.400, 
p/p = (7M? sin? o — 1)/6 
p/p = 6M? sin? o/(5 + M? sin? o) 
T _ (7M? sin? « — 1)(M* sin? o + 5) 
if 36.17? sin? o 


Po _ ( 6M? sin? o ) ( 6 a 
Po M? sin? o + 5 7M? sin? ¢ — 1 


In these equations, p, p, 7, and po denote, respectively, 
free-stream static pressure, density, temperature, and 
stagnation pressure. A dash over a letter indicates 
that the quantity is evaluated immediately behind the 
shock wave. 

The curves W/ sin o = constant are important as 
being level curves for the ratios given above. If this 
family of curves is drawn on the same axes with the 
family of shock polars, the resulting diagram is of con- 
siderable value for the graphical calculation of dynamic, 
as well as kinematic quantities related to oblique shocks. 
The pressure, density, and temperature ratios as func- 
tions of VV sin o for y = 1.400 are shown in Table 1. 
The level values of .V/ sin o are to be preferred, for ex- 
ample, to those of the throttle factor fo/po. The 
simple physical interpretation of \/ sin o is a further 
advantage. 

It seems not to have been observed that the curves 
M sin ¢ = constant are susceptible of a simple para- 
metric representation by means of rational functions. 
The curves of the family may therefore be plotted 
rapidly with an arbitrarily high degree of accuracy. 
For simplicity, the derivation is given for y = 1.400; 
the general case requires only formal changes. 

Received April 26, 1949. 
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CURVES OF CONSTANT .V/ SIN o 
Consider the equation 
Msine =C 
where 


tano = (M* — x)/y (2) 


where x and y are the rectangular coordinates of the 
velocity behind the shock wave and J, taken in the 
x direction, is the quotient of the free-stream velocity 


by the critical velocity. Then,* 


_ (M* — x)(M%x — 1) 


yz = ose (3) 
1 + (5/6)M* — M*x 
This set of equations is completed by the identity 
M? = 5M**/(6 — M*’) (4) 


TABLE 1 


Pressure, Density, and Temperature Ratios as Functions of 





M sin o 
y = 1.400 
M sin o p/p p/p Po/ Po rie 
1.0 1.0000 1.0000 1.000 1.0000 
tA 1.2450 1.1691 0.999 1.0649 
1.2 1.5133 1.3416 0.993 1.1280 
1.3 1.8050 1.5157 0.979 1.1909 
1.4 2.1200 1.6896 0.958 1.2547 
1.5 2.4583 1.8621 0.930 1.3202 
1.6 2.8200 2.0317 0.895 1.3880 
L.7 3.2050 2.1976 0.856 1.4583 
1.8 3.6133 2.3592 0.813 1.5316 
1.9 4.0450 2.5157 0.767 1.6079 
2.0 4.5000 2.6667 0.721 1.6875 
2.1 $.9783 2.8119 0.674 1.7705 
2.2 5.4800 2.9512 0.628 1.8569 
3.3 6.0050 3.0846 0.583 1.9468 
2.4 6.5533 3.2119 0.540 2.0403 
2.5 7.1250 * 3.3333 0.499 2.1375 
2.6 7.7200 3.4490 0.460 2.2383 
2.7 8.3383 3.5590 0.424 2.3429 
2.8 8.9800 3.6636 0.390 2.4512 
2.9 9.6450 3.7629 0.358 2.5632 
3.0 10.3333 3.8571 0.328 2.6790 
3.1 11.0450 3.9466 0.301 2.7986 
3.2 11.7800 4.0315 2.9221 
3.3 12.5383 4.1120 3.0494 
3.4 13.3200 4.1884 3.1803 
3.5 14.1250 4.2609 3.3150 
3.6 14. 953% 4.3296 5- 3.4538 
ee 15.8050 4.3949 G2 3.5963 
3.8 16.6800 4.4568 0.1645 3.7427 
3.9 7.5783 4.5156 0.1510 3.8927 
4.0 18.5000 4.5714 0.1388 4.0469 
4.1 19.4450 4.6245 0.1275 4.2048 
4.2 20.4133 4.6749 0.11738 4.3666 
4.3 21.4050 4.7229 0. L080 4.5322 
4.4 22.4200 4.7685 0.0995 4.7017 
4.5 23.4583 4.8119 0.0917 4.8750 
4.6 24.5200 4.8532 0.0846 5.0524 
4.7 25.6050 4. 8926 0.0781 5. 2334 
4.8 26.7133 4.9301 0.0721 5.4184 
4.9 27.8450 4.9659 0.0697 5.6073 
5.0 29 0000 5.0000 0.0617 5.8000 
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If / and o are eliminated from these equations, there The latter intercept is a maximum—namely, 1—when 


remain the pair of equations 
(C? + 5)M* — 6M*x — 6(C?-— 1) = 0 (5) 


C*y2(6 — M*) = (M* — x)? xX 
[(C? + 5)M* -6C] (6) 


These might be taken as parametric equations of the 
required curves, but one can do better. In Eq. (5), 


place \/*x = t, from which 


M* = 6(t + C? — 1)/(C? + 5) 


and, since x? = f?2/M*, 
x? = (C? + 5)t?/6(t + C? — 1) (7) 


These substitutions are made in Eq. (6), which there- 
upon reduces to 
‘ C? — 1)? (6 — d¢ — 1) 
gras ; 2 — ' (8) 
6C7(¢ + C? — 1) 
Eqs. (7) and (8) constitute the parametric representa- 
While the parameter / 
may be eliminated, the resulting rectangular form is a 
rather unwieldy fourth degree equation. 


tion of the required curves. 


Some obvious properties of the curves may be noted. 
In the construction of a curve, only the positive values 
All the curves are 
All have the x 


of x, of course, are of significance. 
symmetric with respect to the x-axis. 
intercept V 6, which corresponds to t = 6. The other 


x intercept, corresponding tot = 1, is V C2? + 5/CV6. 


V6asC— oa, 
The family of curves fills up a circle. 


C = 1, and it decreases toward 1 
In fact, if 
C— o, the parametric equations become 


x =t/V6, y = (6 — d(t — 1)/6 
The rectangular equation is 


x? + y? — (7x/V6) +1=0 


which will be recognized as the boundary circle of the 
family of shock polars. This is the only curve of the 
family 7 sin ¢ = C which has an axis of symmetry 
parallel to the y-axis. 

The shock polars, and the curves ./ sin ¢ = constant 
are shown in Fig. 1. In engineering applications, it 
has been found desirable, for graphical estimates, to 
have both the family of shock polars and the family 
of curves \/ sin o = constant referred to a system 
of polar coordinates. The velocity magnitude im- 
mediately behind the shock, as well as the deflection 
angle, may thereby be read directly. 
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Inelastic Column Buckling’ 


TUNG-HUA LINt 
Bureau of Aircraft Industry, China 


SUMMARY 


In this paper the buckling of a pin-ended column with small 
initial curvature in the plastic range is investigated. The initial 
curvature is unavoidable in actual columns because of eccentricity 
and nonhomogeneity of materials. A method of calculating the 
deflection of a column with known small initial deflection up to 
the maximum load is shown, and an illustrative example is given. 
From the load-deflection curve, the maximum load is obtained. 
The effect of reversal of strain on the middle part of the column 
on the convex side, when the load and deflection are big as de- 
scribed by Shanley,” ? and the effect of variation of tangent modu- 
lus of elasticity across the section are analyzed. The load at 
which this reversal of strain starts and the method of calculating 
the extent of the reversal of strain on the cross section are shown. 
The theoretical analysis given in this paper shows the different 
factors affecting maximum load and the mechanism of column 
failure. The method of computing the deflection curve at each 
load is shown so that the shortening of the column can be found. 
This may be required for analyzing an indeterminate structure 
with inelastic columns as members. 

The mechanism of the start of deflection of a perfect column 
without initial curvature is explained, and the proposed method 
for calculating the load-deflection curve can also be applied to a 
perfect column and its maximum load can be determined. 


INTRODUCTION 


_ FAILING LOAD of columns has been discussed 
by many authors. Within the elastic range, Euler’s 
equation applies very well. For a column with initial 
curvature, the increase of deflection with load is shown 
in Southwell’s* * well-known method for predicting 
critical load from the load-deflection test of the particu- 
lar column with small initial curvature. This shows 
clearly how the column behaves under various loads in 
the elastic range. 

For the critical load of a perfect column in the in- 
elastic range, two well-known formulas have been pro- 


posed. One is the tangent modulus equation P,,. = 
rE J/L* known as the Engesser theory. The other is 
the double or reduced modulus formula P,,, = 7?E,J/L* 


known as the Considére-Engesser theory or Karman 
theory.> £, is the tangent modulus at that load and 
E, is the reduced modulus shown in reference 5. The 
reduced modulus theory is based on the assumption that 
alter the column reaches the critical uniform stress, it 
bends and the strain on one side increases giving change 
of stress according to the tangent modulus while the 
strain on the other side decreases giving change of 


Received May 9, 1949. 

*The author is greatly indebted to Sir Richard Southwell 
and Prof. J. F. Baker, of Cambridge University, for their most 
kind criticisms and suggestions. 

+ Technical Mission in England, S.M., A.F.R.Ae.S., M.1.A.S. 


stress according to the elastic modulus. The reduced 
modulus formula had been generally considered to be 
the true theoretical solution for perfect columns until 
Shanley raised the point that the column may bend 
simultaneously with increasing axial loading as shown 
in his papers, ““The Column Paradox’’* and ‘Inelastic 
Column Theory.”’! He reached the conclusion that the 
maximum load of a perfect column will be somewhere 
between the tangent modulus load and reduced modulus 
load. There is some range from the tangent modulus 
load to the reduced modulus load, so it is desirable to 
have a method of calculating the maximum load more 
definitely. 

In this paper a complete analysis of an inelastic col- 
umn with known small initial curvature is given and 
an illustrative example is worked out. A method is 
proposed for calculating the deflection curve of an 
inelastic column up to the failing load so as to obtain 
the maximum load and to show how the column behaves 
under different loads. This method can also be applied 
to a perfect column with no initial curvature, so the 
load-deflection curve of a perfect column and its maxi- 
mum load can be obtained. 


THEORETICAL ANALYSIS 
The basic assumptions used for the following analysis 


(1) Plane before loading remains plane 
(2) The fibers of the column act as if 


are as follows: 
after loading. 
they were separate fibers and follow the same stress- 
strain curve as in simple compression and tension tests. 

The symbols used are listed below. The modified 
section referred to in the following is the section derived 
from the geometrical section of the column by multiply- 
ing the width of the geometrical section b by /,//, the 
ratio of the tangent modulus at this point to the elastic 
modulus. The 
width of the modified section }’ equals bE,/E. The 
neutral axis of the modified section passes through its 
Ek, changes back to £ in the region where the 
This will be discussed 


This ratio varies across the section. 


centroid. 
compressive strain decreases. 
in more detail under the section on “‘Reversal of Strain.”’ 


List OF SYMBOLS 


= applied load 
= moment of inertia of column cross section 
I’ = moment of inertia of the modified section about its 
neutral axis 
A = cross-sectional area of the column 
A’ = area of the modified section 


ll 


p radius of gyration of column cross section 
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p’ = radius of gyration of the modified cross section 

h = depth of column section 

hy = distance from the neutral axis of modified section to 
outside fiber on the concave side 

hy = distance from the neutral axis of modified section to 
outside fiber on the convex side 

b = width of column section 

b’ = width of the modified section 

i = column length 

E = Young’s modulus (slope of stress-strain curve in elastic 
range) 

FE, = tangent modulus (local slope of stress-strain curve) 

y = deflection of the column geometrical axis from the line 
of thrust 

yo = initial deflection of the column geometrical axis from 
the line of thrust 

x = distance along the line of thrust measured from one 
end of the column 

6 = deflection of the column at center 

5 = initial deflection of the column at center 

i’) = 5 - 5 

Pe, = 2El,'/L? 

P... = eritical load of the column. It is the load at which the 


column starts to increase deflection without increas 


ing load, ie., dy/dP = oo. Within elastic limit 


Paes ™ Pa 
Mo = moment 
n = distance from the neutral axis of the modified section 
o = stress 
a’ = average stress P/A 
Cea’ @ Pani A 
e€ = strain 
m = distance from the neutral axis of the modified section 


to where de/dP = 0 


Since no actual column is absolutely straight and 
homogeneous, an initial deflection is used to represent 
the initial curvature. For a perfect column, it can be 
equated to zero in the following equations. 

In the elastic range for small deflection, from the 


condition of equilibrium, 


_, (d*v dy 
EI (= - ) 4 Py =0 (1) 
dx? dx? 
Let 
. « NX 
y= z 6, Siu (2) 
. nnx 
vo = p 3 6, Sin (3) 


Substitute Eqs. (2) and (3) in Eq. (1), 


x? «MRK n*n* nTx 
= — 0, Sin - +> 3 > Oo» sin +- 
I I L? 


R } y Vt 
ag . NTXx 
— > 3, mn = 0 
EI L 
, . nnx : : 
Multiply by sin dx and integrate from x = 0 to 
x = L; all the terms, except m = n, vanish. 
. nx? - nn? _ ¥ 
—dn + 6, - = + 6,— =0 (4) 
L? L? El 
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6, (; - a 43, n> =0 
EI L? L? 





- 5b» 
6, = F 
1—|[P/(n?r2kI)/L? 
Po = El / iL 
A bn 
On h 
| — (P/n?P.,.) 
cai . = ( l ) s 
0 co. ~ 0. = — | 0, 
l — (P/n*P., ) 
> c 
é./ 4 6,f = 0, 
n*P., — P (n*P.,/P) — 1 
6,’ 61 7 
l wt i 
(P.,/P) — 1 


As P approaches P,,., 6; is much more amplified than 
the amplitudes of other harmonics. 


m =“ ‘ 01 
o’ = pe By = 6,’ = § 


This shows the relation of the increase of deflection 
with load in the elastic range. 

When the column is loaded beyond the elastic limit, 
the flexural rigidities of the different sections of the 
column vary with the load. At a particular section of 
the column under a given load and moment, from the 


condition of equilibrium, 


> 


r= fou q 


« 


* 


M = fon dA (10 


For a small increment of load AP and a small incre 
ment of moment AJ/, from Eqs. (9) and (10), 


AP = f dA (11 


AM = fam dA (12) 


Ac = E,Ae (13 
From the assumption that planes before loading re- 
main planes after loading, 


e= e+ Kn (14) 


where A is a constant. 

Assume the value of £, as constant in the interval 
from P to P + AP and from M to M+ AM. The 
smaller the value of AP and AM, the more accurate 
will be the result. With /, remaining the same, there 
is no change in the modified section within the interval. 


Ae = Ae + nAK (15) 


From Eqs. (11), (15), and (12) 
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INELASTIC COL 
p= fi £,(Aeo + nAK )bdn (16) 
AM = fe (Aey + nAK )bndn (17) 


By definition of the neutral axis of the modified section 


Af en = 0 (18) 


I: ,bdn = EA’ 19) 


f ana = Py’ (20) 


From Eqs. (16) to (20), 
AP = EA’ he, 21) 
AM = EIl'AK (22) 
From Eqs. (15), (21), and (22) 
pa AP 4 AM (2 
EA’ EI’ 


When the increment of load and moment approach 
zero, the increments Ae, AP, and A.JV/ are replaced by 
differentials de, dP, and dA. 


dP dM 


le = + (24 

wodnee Youle | 

k#dP EkdM - 

do = n (25) 
EA’ Fm fed 


Let the column load P be increased to P + AP and 
the deflection y be increased to y + Ay. From the 


condition of equilibrium, A.J = (P + AP)(y + Ay) - 


Py 
AM = —EI'A(d*y/dx?) = (P + AP)Ay + APy (26) 
ay > > 2 
sO Fe ei lk 
dx? iI’ P pr’ 
d*Ay 


2 ] > 
P+aP 
dx? rl oa 


Effect of Eccentricity 


Due to the fact that with deflection the strain, and 
hence the tangent modulus across the section, varies 
from the concave to the convex side, the neutral axis 
of the modified section does not coincide with the 
geometrical axis and changes with load. But the incre- 
ment of load still applies at the geometrical axis and 
hence will have an eccentricity x» on the modified sec- 
tion. This eccentricity increases rapidly as reversal 
of strain occurs. So in the last term in Eqs. (26) to 
(28), y should be measured from the neutral axis of the 
modified section. But y was taken as the perpendicular 
distance from the line of thrust to the geometrical axis, 
so an eccentricity x9 should be added to y in the last 
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term in Eqs. (26) and (28) and these equations become 


AM = (P + AP)Ay + AP(y + x) = El’AK (29) 
12 , P 2 2 
—_ +” a+ & (y+tx.) =0 (30) 
dx* EI’ EI’ 
Integrate Eq. (30) and let the lower limit be x = L/2: 
, "x Dp ) 
dAy — / P+ AF Ay de = 
dx stn i 
7. Peal 
(y+ x)dx+C (31) 
ae 3 hid 
Atx = L/2, dAy/dx = 0, so C = 0. 
2 | 2 
-f fa : - Ay dx dx — 
we 
J. ie ‘3 vo )\dx dx + Ci 32) 
L/2 fa" 
At x = L/2, Ay = Aé, so CG, = Aé. 
“x Dp P 
Aé - f 3 : Ay dx dx 
I 
(y+ x)dx dx (33 


f r AP 
2J L/2 EI’ 
The above integrals are the moment of the area 
under the curves of (P + AP)Ay, EJ’ and AP(y + x 

EI’ between x = L/2 x about the section 
x = x. The Ay curve can be calculated by 
Aé and finding Ay by starting from the center section 
x =L 2. If the Ay at x = Lequals zero the assumed 
Aé is correct, otherwise Aé should be modified until 
this condition is satisfied. For further increment of 
load, the new EI’ at different sections should be calcu- 


lated and the above process repeated. 


and x = 
assuming 


Reversal of Strain 


A typical compressive stress-strain curve is shown 
When the strain is increased from zero, the 
If the strain is increased 


in Fig. 1. 
stress follows the curve OCAB. 
up to A and then reduced, the stress-strain curve will 


not follow A to C, but A to D which has the same 
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Stress-strain curve for mild steel for composition (from 
reference 7). 


Fic. 2 


So be- 
yond the proportional limit, if the compressive strain 


slope giving Young’s modulus of elasticity /. 


of a portion of the section decreases with the increase 
of load, & would occur instead of /;. So where de/dP 
> 0, E, would occur, and where de/dP < 0, / would 
occur. 

Now let us consider the distribution of strain across 
the section when the stress is beyond the proportional 
limit. 


AP AM 


n (24) 
EI’ 


Ae = _ 
FA’ 


AM = (P + AP)Ay + (vy + x)AP (29) 


From the above two equations 


AP , (P+ AP)Ay + (y + x0)AP | 
a T 


Ae = - 34) 
EA’ EI’ 

ae _ sit P(Ay/AP) + y + Ay + xo » (35) 

AP EA’ EI’ 


When Ae/AP = 0, 
] P(Ay/AP) + y + Ay + Xo 
= — —— y 


A’ r 


11 


SCIENCES—MARCH, 1950 


ne : (36 
P(Ay/AP) + y + Ay + Xo 

As the value of AP approaches zero, Ay approaches 

zero and Ay/AP = dy dP, then Eq. (36) becomes 

2 
= 2 (37 
P(dy/dP) + y + Xo 

m is the point on the convex side, beyond which th: 
compressive strain decreases and /, changes back to £ 
When 7; is greater than /», there is no region of reversal 
of strain. Up to m = /y, the tangent modulus is effec 
tive across the whole section. Further load and deflec. 
tion would introduce a double modulus, by which js 
meant that Young’s modulus F extends over a portion 
of the section from the convex side and the tangent 
modulus covers the remaining portion. As load and 
deflection increase, there will be an increase of the region 
where Young’s modulus is valid. The modified section 
at different stages of loading is illustrated in Fig. 3, 
which is based on the linear variation of /, with e that 
occurs in a region beyond the proportional limit (as 
shown in Fig. 2) for most engineering materials. For 
regions in which this Inearity does not hold, the modi- 
fied section should be based on the actual variation of 
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| j Se. 86. 
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; FiG 3e 
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= 31h 


Modified sections of mid-section at different loads. 
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AT o's 2900 Ky /om? Substitute the above in Eq. (50) 
(56 — i 
—_ n » » 
- FIG 44 ee sc P+AP _  . ATX. 
ae 1: & cou Aé, sin + >°A6, sin + 
f L? I EI’ I 
/PTOaches L, = .0439h AJA = .82 a 3 : . 
me . ' AP nTx nTx 
1S 1! ih " ° 7. ° we 
. 2 f= p>—_2% (Xs, sin + Sea, sin ) = 0 
‘yr — 
ne 2 =z - 30h kil L L 
(27 } 
\o Son ‘oC . . ; 
—| I’ at different sections is assumed equal to /;’ at 
hict center section throughout the following analysis. 
1ich the Fig 4b , : , 3 
k th: - a Multiply the above equation by sin (mmx L)dx and 
LCK to F ~ ° ° — > 
E a X,=.0354h A/a = .82k integrate from the limit of x = 0 tox = L. All the 
reversal i | ; ' ae ; = ci 
is eff y 1/1 = .842 2, 2 leh terms m # n vanish and the terms m = m give the 
1S effec. 7 ° e . 
” h,- 7. = + 30h following equations: 
d deflec. . 
” Suh 9 9 
which jg =r nts? P+AP.. a .. 
os — Aé, + 7" Aé, + — - (6, +a,) = 0 
portion Fig 46 L? EI, EI, 
tangent S35 : 
: | t : n*n? EI’ . " . —_ 
vad and 2b 1376 X,- -028h - —_ + (P + AP) |A6, + AP(6, + a,) = 0 
€ region i + AJA =-796 _1'/I = 524 L* 
sant ‘ - 3S} 
Section _ 47. : “ (= — , 
Fig. 3, Substitute P.,, for r?/I,’/ L?, 


1 é that 
mit (as FiG 4d (n*P.,, — P — AP)Ab, = AP(6, + a,) 


| -20 h 
Ss. For L . =~ % Aé 6, +a 
> modi- Teo — <= - = a - (39 
t a AJA: .% 1/1 = .75 AP n’P.,, — P — AP 


tion of 








46 _ > Ab, _ 6, +a, 
AP AP n*.,, — P — AP 
, F6.de As in the elastic case, the increase of first harmonic is 
x-+O 


egies much more than the other harmonics as P approaches 
- 74 b = 


AYA =.74G Wi ».146 P..,.;, so the increments of other harmonics are neglected. 











Aé >" 6b» + Qn a, 6; + ay 


Fic. 4. Modified sections at different sections of column = ~ 


AP “n?P.,—-P—AP n°P.,—P— AP 


; = c ; As stated above, the J’ at different sections are 
E, with e. The reversal of strain would affect both the ’ ae — 
assumed to be equal to /,;’ at center section. It is fur- 


! value of /’ and the position of neutral axis. ; : s ; 
# ; : pe ; ‘ ; ther assumed that the amplitude of the first harmonic 
With given form of section, stress-strain curve of the 
material and the known initial deflection (the initial 


deflection may be found from Southwell’s* * test in 








elastic range), the load-deflection curve may be obtained ] 

by Eq. (7) until the extreme fiber stress on the concave me eS ie a aa aS eS 

side is equal to the elastic limit of the material. Then | |e 

for subsequent increase of load, get the first approximate l 4] | =. | | a 4 
curve of the increment of deflection Ay as if the column 2 ee ae ee ne ae eee eee 


were of constant //’ equal to that at center section. 


From Eq. (30) ae a el sain 











d*Ay (P + AP) AP 8 Lo 
: - sen Ay + — (y+ %) = 0 (30) 

dx? EI El Oth 
Let 02 | 
a MEX 2. Le, OS ais 

dy = > 46, sin-——, y= > 5, sin — " 3 

L i 
i 4 4 O4 
4 > s nNTX — 
s ad, Sin I ¥% ar at': 2900 kg/cm? : 
; + ’ P "OS + 
d*Ay ne? g NHX Ps Sos 
~ aie id p ~ Ad, sin Fic. 5. Deflection and eccentricity of modified sections along 


1s. dx? Be L length of column 
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| -“ 
— een -O03h 
Ba 
= = sue — ay 
A 
+>} yx 002 h 
FIRST APPROXIMATE Curve. 4 Y_= 00336 sin oh 
i x CALCULATED vA.ueS BY Eq. 33. 
a = V4 a ee ee = ime 
° a ome 3 4 5 .G 7 8 9 1.0 
| x-£ 
Fic. 6. The approximate and calculated curve of the increment of deflection due to 40’ = 10 kg. per sq.cm. at o’ = 2,900 kg 
per sq.cm. 


of eccentricity a, and that of the deflection curve 6, are 
equal to their values x» and 6 at the center, respectively. 


Aé6 _ 5 + Xo (40) 
AP Pa, —P— AP 
As AP approaches zero, A6/AP = dé/dP. 
dé = (6 + Xo) _ (41) 
dP 0 ee 
From Eq. (37), atx = L/2 
ae. en 
P(di/dP) + 6 + xo 
From the above equation and Eq. (41) 
a. an 
(6 + x0)[(P/Pe, — P) + 1] 
ur (42) 


7 (6 + Xo) [(Per., Pax aa P)| 


From Eq. (42) and the relation P,,, = 7*EI,’/L?, the 
modified section at x = L/2 for the first approximation 
can be determined and so the values of A,’, J)’, and 
dé/dP at the center section are obtained. With the 
above assumptions, the curve of the increment of de- 
flection for the first approximation is first harmonic. 


Ay = Aé sin (rx/L) 
Ay/AP = (A6/AP) sin (2x/L) 
dy/dP = (dé/dP) sin (rx/L) 


With this dy/dP value put in Eq. (37), m is found for 
different sections. If m, > A» there is no reversal of 
strain. If n, < hy reversal of strain occurs in the 
region n > 7, ton = fe. If the deviation of J’ at various 
sections from /,’ is negligible, no revision of the incre- 


ment deflection is necessary. If the deviation is not 
negligible, apply Eq. (33) and by the moment area 
method, obtain the increment deflection curve due to 
the increased bending moment. If this curve deviates 
negligibly from the first approximate curve, the first 
approximate curve is considered to be accurate enough. 
If not, the increment deflection curve obtained from 
the bending moment is used as an assumed second ap- 
proximation and the process is repeated until the curve 
obtained from the increased bending moment agrees 
closely with the assumed curve. After the Ay is ob- 
tained in this way, the increment of strain across the 
section can be obtained by finding Ae) from Eq. (21). 
Where there is no reversal of strain across the section 
and the £, vs. e curve is linear, é) can be taken from the 
stress-strain curve corresponding to the average stress 
of the section P/A. From Eq. (29) 


a « (P + AP)Ay + (vy + xo) AP 
BY’ 
Then from the strain 


where AK equals A(e; — e2)/h. 
Then add 


distribution, the new //’ can be computed. 
the further increment of load, and repeat the process 
until the increase of deflection becomes very large indi- 
cating that the buckling load is being approached. 
This method involves the assumption that in the inte!- 
val from P to P + AP, EI’ remains constant at the 
The error in this assump- 
AP taken should be 
The con- 


value corresponding to P. 
tion diminishes as AP decreases. 
varies sharply with strain e. 


1 


small where F, 
vergence of the above successive approximation method 
is shown in Appendix (II). 

For an actual column, the value of E, decreases al- 
most linearly with strain beyond the proportional limit 
as shown in Fig. 2. The value of J’ actually varies little 
from the ends to the center of the column before the 
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reversal of strain occurs. For a first approximation, 
assume EI’ to be constant and equal to EI,’ at the 
center of the column. For a steel column of rectangular 
section [calculated in Appendix (I)], the difference be- 
tween the value of E,J at the ends and £/,’ at the 
center is about 1 per cent before the occurrence of the 
reversal of strain, so the above first approximation curve 
is close enough for the calculation of the increment of 
deflection until the reversal of strain occurs. After the 
reversal of strain occurs, the increment of deflection is 
calculated by the moment area method based on Eq. 
(33). For those sections with reversal of strain, Eq. 
(21) should be used for the increment of axial strain ép. 

The load-deflection curve of a steel column of rec- 
tangular section, of slenderness ratio of 75 and with 
initial deflection of 0.001 in Karman’s test,® is calcu- 
lated by the above method and shown in Appendix (1). 
The calculated curve agrees very well with tested 
values. 

For a column of small L/p, as load increases, E/ re- 
mains constant up to the proportional limit. Beyond 
the proportional limit, // decreases due to the reduction 
of EI to El (2, < #). Then with increase of deflec- 
tion, a portion of the cross section has negative de/dP 
at the convex side and double modulus occurs. For the 
portion 7 > m, @, changes back to £, but on the rest 
of the section, F, still exists and decreases according 
to the stress-strain curve of the material. Generally 
as the reversal of strain occurs, the effective moment of 
This increase of ec- 
Even- 


inertia and eccentricity increase. 
centricity increases deflection as load increases. 
tually the deflection increases indefinitely for a definite 
increase of load and the maximum load is reached. 

A perfect column will remain straight up to the 
tangent modulus load. If there is no slight disturbance, 
the column may remain straight even beyond the tan- 
gent modulus load. A slight disturbance is assumed to 
be present all the time. This slight disturbance is here 
represented by an infinitesimal deflection 6 = « = 0. 
When this 6 = ¢ is very small, EJ’ may be taken as 
constant along the column equal to £,/. From Eqs. 
(41) and (42) 


d6/dP = (6 + Xo/P.e., — P) (41) 


9 


m = a (42) 
(6 + X0)(P. /Po. — P) 
até = <0, x% = 0 
di/dP = (¢/P.,. — P) 
nm = —_- 
= : 
dP s/t oe. — £) 


Here P,,, = 7°E,J/L*. When P,,. — P approaches the 
magnitude of e, d6/dP increases to a finite value. Let 
P.,. — P = ne, where n is a finite number. 


di/dP = 1/n 
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—np’ * —np’ * 
m 5 ai . , 49 ° 
P.. (x2EA 'p’2/L?) 
m = — nL?/r?EA’ (43) 


Equate 7 to —/p. 
ho = nL?/r?EA’ 


Py. — P _ wEA'hy _ 


i ae a _ T*E Aho (44) 
€ I? L? 

As P approaches P,,., at first, m is big and m > he 
until the ratio (P,,. — P)/e reaches the value shown in 
Eq. (44); then 7, = fz. Further decrease of m occurs as 
P approaches nearer to P,,. and will be less than hp 
and double modulus occurs. Where reversal of strain 
occurs, the tangent modulus is changed back to elastic 
modulus, so /’ and P,,. increase im value and the eccen- 
tricity of the modified section increases rapidly to finite 
value. Then 6 increases with additional load according 
to Eq. (41), mainly because of x». When 6 increases to 
a finite value, the variation of tangent modulus across 
the section will also increase eccentricity of the modified 
section. 7», decreases according to Eq. (42) and the 
eccentricity of the modified section increases. So the 
deflection increases with further additional load. With 
increase of deflection, the compressive strain on the 
concave side increases more than at the neutral axis of 
the modified section. The tangent modulus on the 
concave side decreases more than that at the neutral 
axis. (For all engineering materials the tangent modu- 
lus decreases with strain beyond elastic limit.) This 
decrease of tangent modulus will have the effect of 
reducing the moment of inertia of the modified section 
and hence the critical load. The combined effect of the 
reversal of strain, the decrease of tangent modulus at 
the neutral axis and the further decrease of the tangent 
modulus on the concave side will give a critical load at 
sach value of load. Deflection starts at no reversal of 
strain 7; > fe and not at m, = O as in the reduced 
modulus theory. This explains why deflection starts 
at the tangent modulus load and not at the reduced 
modulus load for a perfect column with no initial curva- 
ture. With further increment of load and deflection, 
n decreases and approaches zero at critical load, while 
the compressive strain increases more on the concave 
side than at the neutral axis. Since the tangent modulus 
decreases with compressive strain beyond the propor- 
tional limit, the decrease of tangent modulus is more 
on the concave side than at the neutral axis of the modi- 
The EI’ of the modified section is thus 
Hence 


fied section. 
less than E,J in the reduced modulus theory. 
the reduced modulus load is not reached. Taking into 
consideration the effect that n; decreases gradually from 
h, and P,, and x» increase with the decrease of m, the 
increment of deflection with load beyond the tangent 
modulus load can be found in the same way as for the 
column with initial curvature, as explained previously. 
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Fic. 7. Load-deflection curve for steel column L/P = 75. 


Curve from test in reference 7. 
deflection § 0.001h. 


x calculated value. Initial 


With load-deflection curve computed, the maximum 
load can be obtained quite closely. 

The critical load would depend on £/’ along the 
whole length of the column. For an actual column with 
small initial curvature, when the load increases near 
the critical load, the deflection increases and the de- 
nominator of the right side of Eq. (37), especially the 
term dy/dP, increases and » decreases, approaching 
zero at critical load. This decrease of 7; increases //’ 
of the modified sections. On the other hand, with in- 
creasing deflection, the compressive strain increases 
more on the concave side than at the neutral axis. 
So the tangent modulus decreases with the distance 
from the neutral axis on the concave side. This causes 
the tangent modulus on the concave side to be less 
than the tangent modulus at tangent modulus load. 
The excess of the effect caused by reversal of strain over 
the effect caused by the decrease of tangent modulus 
on the concave side at the maximum load, increases 
the maximum load above the tangent modulus load. 


50 
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Fic. 8. Column theories and test data (Van Den Brook), 


reference 1. 
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Ato’ = 2,400, 6’/h = 0.0016. 


SCIENCES—MARCH, 1950 

For very short columns of 24ST, the rate of the decreag, 
of tangent modulus with increase of strain near th, 
critical load is much less, so the effect of the decrease 
of tangent modulus on the concave side on //’ is much 
‘ess. The tangent modulus at the neutral axis is quite 
near to the tangent modulus on the concave side in the 
reduced modulus theory. And in very short columns, 
the tangent modulus does not vary much with the 
distance from the neutral axis and would be quite close 
to the value at neutral axis and to the tangent modulys 
value of the reduced modulus theory. The difference 
in value of the tangent modulus and elastic modulus 
increases with the critical stress, so for short columns 
with high critical stress, the effect of reversal of strain 
is more prominent. Hence, the maximum load would 
deviate more from the tangent modulus load. Test 
points of the very short columns, therefore, often lie 
closer to the reduced than to the tangent modulus 
curve as shown in Fig. 8 and other similar tests. 
APPENDIX (I)—CALCULATION OF LOAD-DEFLECTION 
CURVE OF A STEEL COLUMN 


In the following, the load-deflection curve of a steel 
column tested by Karman,° of rectangular section with 
initial deflection of 0.0014 at the center and slenderness 
ratio L/p of 75, is calculated based on the method 
shown. The following data are taken from reference 5. 
The stress-strain curve of the material is taken from 


Fig. 6 of the reference. The tangent modulus is calcu- 


lated from the stress-strain curve and plotted in 
Fig. 2. 
Lis = ta, FE = 2,170,000 kg. per sq. cm. 
5 = 0.001h 
At o’ = 34,000 Ibs. per sq.in. or 2,400 kg. per sq. cm., 


6’ = 0.0016h. For this deflection and load, the column 
is in the elastic range. So Eq. (8) is used to find 6). 


6’ = Or 
(P.,./P)—-1 
P.. = Ae...’ P.= Ao’ 

5’ = 01 (45 

(Ger.’/o’) — 1 
; nr El rE ° 
Co, = = (46) 

L°?A (L/p)? 


“ts 2,170,000 
= ae = 3,810 kg. per sq.cm. 


io” 


From Eq. (45) 


6 
(3,810/2,400) — | 


0.0016h 


5,/h = 0.000938 


The first harmonic is much more amplified than other 


harmonics, so the deflected curve, at o’ = 2,400 kg. 
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TABLE 2 


Calculation of Section Properties of Modified Sections—at o’ = 2,900 Kg. per Sq.Cm., at Different Values of x 
Based on the First Approximate Curve Ay = Aé sin (rx/L) 


«/L 1/2 1/3 1/6 1/10 1/20 (0 


o’ Kg. per sq.cm. 2,900 2,900 2,900 2,900 2,900 2,900 
eo Fig. 2 0.00142 0.00142 0.00142 0.00142 0.00142 0.00142 
y'/h 0.01024 sin (rx/L) 0.01024 0.00888 0.00512 0.00317 0.0025 0 
nh? y’ " ’ ee ei 
el — & 0.000215 0.000186 0.000107 0. 000066 0.000052 0) 
[* hk 
a — & e m ” } 5 
i“ oii 1 h; 0.00153 0.00151 0.00147 0.00146 0.00145 0.00142 
U 
Qa — ee = 
C2 eo + j — hy 0.00182 0.001388 0.00137 0.001389 0.00140 0.00142 
1 
E,,/10® Kg. per sq.cm., from Fig. 2 1.40 1.45 1.55 1.55 1.60 1.62 
Ey,/10®8 Kg. per sq.cm., from Fig. 2 1.85 1.80 1.75 1.70 1.67 1.62 
E/E 0.645 0.668 0.714 0.718 0.737 0.746 
E,2/E . 0.783 0.770 0.746 
Xo/h Fig. 4 0.0465 0.04389 0.0354 0.0282 0.008 0 
A’'/A Fig. 4 0.81 0.82 0.821 0.796 0.750 0.746 
I'/bh3 ‘Fig. 4 0.0674 0.0685 0.0702 0. 0686 0.0625 (0.0622 
he — m 
———~ Fig, 4 0.37 0.36 0.30 0.20 
h 
Gy = C2 
en pee =m 0.00140 0.00140 0.00141 0.00141 
1 
Etm,/10® Fig. 2 1.68 1.67 1.675 1.675 
E 
7 0.774 0.770 0.772 0.772 
E 
m/h Fig. 4 0.0835 0.096 0.165 
dy = = 
> —= 0.941 0.815 0.471 0.291 0.147 
dP 
y/ht 0.01124 0.00975 0.00562 0.00348 0.00175 0 
Modified section Fig. 3e Fig. 4a Fig. 4b Fig. 4c Fig. 4d Fig. 4e 
dé 5 + Xo 0.01124 + 0.0465 ; 
a = P = — = 0.941; P(dy/dP) = 0.941 sin (ax/L) 
dP Poe = 2 (3,078 — 2,900) /2,900 . 
ty = é6sin (wx/L) = 0.01124 sin (rx/L). 
per sq.cm. would be very close to the first harmonic*™™ (€: — @2)/h = —6x?/L? 
curve. ev ssrs 
" A(e; — €2) k= —A6é'r? L? (48) 
5, tee. d*y Tr TX er : ae 
y = 6sin —, -_=-— 6 sin Dividing the denominators of Eq. (40) by A on both 
1a dx . L? i, sides, 
; + 4) ° n 
From Eq. (14), Ab _ 5+ x “ 
e= e+ Kn, €; — & = Kh, A(e; — @) = hAK Ao’ o.,,. — a’ — Ao’ 
Then consider a segment of beam of uniform curva- From Eq. (42), with P.,, = Age, / and P = Ao’, 
ture with radius R and extending an angle #. Equate 
the difference in length of the extreme fibers of this —p’? (30) 
° ° ° i = - ie ’ 
segment to the difference in strain at the extreme fibers (6 + x0) (Ger.,’/Ger.,’ — 0’) 


times the length of the neutral axis of the segment. 
Eqs. (49) and (50) are based on the assumption that 


(R + )@ — (R — he)@ = RO(e, — e2) EI’ is constant throughout the length of the column, so 
h = Rie; — e) these equations are only justified when the variation of 

: EI’ is shown to be negligible. Eq. (50) is used for first 

é: — é: = h/R = Kh approximate calculation of m in Table 1, row 23, but 

K = 1/R = d*v/dx? this value is to be revised if later successive approxima- 


, , THe 7 tion shows necessary revision. 
(€; — @2)/h = d*y/dx* = —6(2?/L?) sin (rx/L) (47) : : , P 
, . Eq. (49) is used to calculate the deflection curve ol 


Atx = L/2 the steel column until the reversal of strain occurs. 
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Before the occurrence of the reversal of strain, the value 
of EI’ along the length deviates very little from its 
value at the center section. This is shown in the close- 
ness of the ratio HJ EI’ to unity as listed in Table 1, 
row 15. 

So calculation of the deflection curve up to a’ = 2,900 
kg. per sq.cm. in Table 1 based on the assumption of 
constant //’ along the length of the column at each load 
is justified, and the deflection curve at ¢’ = 2,900 kg. 
per sq.cm. is very near to first harmonic. After the 
occurrence of the reversal of strain at o’ = 2,900 kg. 
per sq.cm., the ratio El EI’ is 0.912. Tts deviation 
from unity would suggest the use of moment area 
method by considering the variation of //’ along the 
length as shown in Tables 2 and 3. 

As mentioned above, the deflection curve at o’ = 
2,900 kg. per sq.cm. is very near to first harmonic, so 


y = 6 sin (rxv/L) 
and 


(d*y/dx?) = — (r?/L?)6 sin (rx/L) 


From Eq. (48), A(e: — e2) 4 and then (e; — e)/h for 
different sections can be found. With é obtained from 
the stress-strain curve corresponding to o’, the values 
of e; and e: can be obtained by Eq. (14) with assumed 
eccentricity and later checked with the eccentricity 
found from the modified section. Then /,, and /;, are 
read from Fig. 2 corresponding to e; and és. 

In order to compute n;, dy/dP is required. So, for 
the first approximation, assume the column to be of 
constant EJ’ equal to that at the center section and 
assume the first harmonic amplitude of the (vy + x») 
curve at this load equal to its value at center section 
(6 + x»). Assume ny; and draw the modified section. 
Check the value of » with Eq. (50). Revise m, if it 
does not agree with Eq. (50). 

With /,’ and p’ ° known, calculate o,,.,’ by the equa- 
tion o,,,' = mI’ L?A. Apply Eq. (49) to find Aé. 
The increment deflection curve, based on the above 
assumptions, is a sine curve, Ay = Aé sin (rx /L), so 
dy dP can be readily found by multiplying the value of 
dé,dP by sin (rx L). Withdy dP inserted in Eq. (37), 
the value of n, for other sections where reversal of strain 
exists can be found in the same way. The values of /’ 
and x, at different sections are computed in Table 2 
and plotted in Fig. 5. 

With arbitrary Ao’, the first approximate value of 
Aé can be obtained from Eq. (49) and the first approxi- 
mate increment deflection curve would be Ay = Aé 
sin (7x/L). Use this curve as first approximate Ay curve, 
calculate the deflection at various sections because of 
the increment of bending moment PAy + AP(y + Xo) 
as shown in Table 3. If the new curve does not agree 
with the assumed Ay curve, the new curve should be 
used and dy/dP obtained by assuming it equal to Ay/AP 


from the new curve and 7 recomputed and x», J’ t¢. 
vised. Then compute the deflection curve caused py 
the increment of bending moment PAy + AP(y + y, 
as before. This increment deflection curve would be 
the second approximation. But as shown in Fig, ¢ 
the computed curve agrees extremely well with the 
assumed sine curve, so no further step is necessary. 
From Fig. 7, the curve is quite flat at this load, so this 
load may be taken as the maximum load. In general, 
this method of successive approximation converges 
rapidly to the final Ay curve. The convergence of this 
method is shown in Appendix (II). 


APPENDIX (IT) 


The above method of successive approximation js 
shown to be convergent as follows. From Eq. (33), 


“a a > > 
Ay = Aé — / | rts Ay dx dx — 
L/2J L/2 El’ : 
‘x 7x Ap 
J | . (vy + Xo)dx dx 
L/2J) L/2 EI’ 


as > tat 
‘im | J ! (PP + AP)Ay + 
L/2J) L/2 EI’ 


AP(y + xo)|dx dx (51 


is constant. P and y, which is a function of x, are 
known. The column is divided into n equal segments 
and they are numbered from the center of the column 
from | to m on each side. Let f, F, G, ¥, and @ denote 
functions. AP is given. From Eq. (33), 


Ay; = fi(A6, Ii’, Xo,) 


Aye fo( M6, [,', J,’ Xo,» Xo.) 


II 


Ay; = fi(A6, i’, I,’, ce ey Be XO,» XO. «~~» XO) 
T,’ and xo, are functions of Ay,. 
I,’ = y(Ayi) 


Xo; = i(Ay;) 


Ay. = filAd6, Yi(Ay1), d:(Ay1)] = Gi(Ad) 

Aye = fo{Ad, Yi[Gi(Ad)], Yeo(Aye), $1[Gi(Ad)], (Aye 
= G»(Aéd) 

Hence 


Ay, = G;(A8) 


dl,’ dl,’ dAy; - a7,’ 


dAé dAy; dAé dAé 


From Eq. (51) 


Aé = F(i,', De... 7 PP Ay, Ayo, secs Ayn; X09 VO.» om 
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INELASTIC 


dAé = 
x( OF dl,’ 


OF dAy, ; OF we) das (52) 
i\O/;’ dAé 


OAy; dA6d OX; dAéb 


If the above condition is satisfied, increment of deflec 
tion may increase without the increase of AP. Since 
Pp + AP is less than the critical load, the moment due 
to the increase of deflection alone is not sufficient to 


produce that deflection. Therefore 


oF day, 


—_ oe 
OAy, ddd ; 


: (= dl, 3) 
Oxo; dAd ‘ 


> 
i-1 \Ol,’ dAé 
In the parentheses, the second and the third terms are 
positive because of the fact that the moment area in- 
creases as the deflection and eccentricity increase and 
that the eccentricity increases with deflection and 
day/dP. But the first term may be negative when /’ 
increases at the sections where reversal of strain occurs. 
Considering the extreme case with reversal of strain 
throughout the length of the column, from Eqs. (36) 


and (51), 


~ fs . 
Aé = / | - | P + AP) a + (y+ | dx dx 
J it2J L/2 EI’ AP 


9 
s2 


p 


nailes (P + AP)(Ay/AP) + y + xo 


I' = A'p’* = —A'm[(P + AP)(Ay/AP) + vy + xo] 


re *x >) 
Aé = | | ad dx dx (54) 
L/2 JL/2 —EA'n 
dAé = (= = a-A ™) dAéb (55) 


10(—A’m); dad 
Let z in the following be the distance from the ex- 
treme fiber on the convex side, 2, = fe — (—m), m = 
E/E, and bE,’/E = f(z), where E,’ 
modulus without considering the reversal of strain. 


Let 
Q= / of(z)dz 
J 0 . 


and Q’ be the moment of the modified section area 


is the tangent 


about z = 0. 


+ > 


A’ = / mfteide +f f(z)dz 
J 0 2 


7 7) 


21 ha 
Q’ = / mf(s)2dz + / f(z)s dz 
0 4 


« - 


he = 4 + (—m) = Q'/A’ 


*2) mh 
= / (m — 1)af(z)dz + / 2f(z)dz — 
0 0 
72) *h 
21 fi (m — 1)f(s)\dz — x J f(z)dz 
0 0 
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—A’'n = 21 / (m — 1)f(z)dz — / / (m—1)X 
J 0 J0 JS0 
f(z)dzdz + Q — af (m — I1)f(z)dz 
0 


ba / S(z)dz 
J 0 
Q - af f(s)\dz — 
Jo 
/ / (m — I)f(s)\dzdz (56) 
0 J? 


From the above equation, |d(—A’m)/dz,| < 0, reversal 
of strain increases with the increment of deflection, so 
(dA6/dz,) > 0. 


d(—A’‘m) d(A’n) dz, os 
= < 0 (4) 
dAd dz, dAé 
From Eq. (54) 
[0A6/0(—A’m);| < 0 (58) 
OA6 d(—A’n); ‘ 
> 0 (99) 
0(—A’n), dad 
— ons dA’ 
a Se (60) 


i=10(A ‘m) i dAé 


Equate the right side of Eq. (52) to that of Eq. (55). 


From Eqs. (53) and (60) 


po tlt OF day oF &s) 
0< > : <1 
i=1 \O/;’ dAé OAy; dAd Ox»; dA6 


(61) 


Let the subscript I, II, III, IV, , NV denote the 


order of approximation. In the successive approxima- 


tion method shown in this paper, the first increment de- 
flection curve is obtained by assuming it to be first 
harmonic, 

Ay; = Ad, sin (rx /L) 
Then from Ay,/AP, get m, Xx» and then with Ay, curve, 
calculate A.\//E/’ with the new 7; and x». Then 


af "3 
nl 
Ady; = / / . f dx dx 
J L/2) L/2 EI’ 
= 4* 
Ayy =A 64; — / / — dx dx 
F L/2J) L/2 EI’ 


3ased on Ayyy, calculate 
This 


Then from Ay,,/AP, get m, Xo. 
Ayr, and Ady; by the moment area method. 
process is repeated until Ad, is close enough to Ad, 


The above can be expressed as follows: 


Aéry = Ad; a >| 


i l 


OFT 
ol,’ 


li’ —TJ,’)+ 


oF =) ” 
(« , = Xe, SN 62 
OX L 


where xo, is the eccentricity at the center section in the 
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first approximation. 
curve based on Ay,; is used, while in computing Aé;; 
the moment curve based on Ay, is used. 

‘(OF dl,’ . OF dxo, 
Adi = Ady + Zz Fi = 

t Ol; dAé6 OX0; dAé 


=1 


‘—" (OF dAy; 
(Ady; — Ad;) + 2) (as — Aéd;) 
. , » im dAé “ 
OF ddy; 


oo (= ai,’ 
OAy; dAé 


Abi, — Ad, = > 
i - i=1 \O/;,’ dAé 
a) (Aé;; — Ad) (63) 
OX; dAé 
Here Aé is the final value to be found. 


As = 05, + ¥ | Seay —I,')+ 
i=1|.O7;,' 


1 


oF (2 — Xo, sin *) + 
OX: ‘6 


‘=! OF day, 


(A6 — Aéd;) (64) 
i=1 OAy; dA6é 
From Eqs. (62) and (64) 
t n " % F I, 
Ab — Ab, = (= di OF ddy re 
i=1 \O/,’ dAé OAy; dAé 
wl ou) (AS — Aé;) (65) 
OX0; dAé 
From Eqs. (63) and (65), 
Aé as Adérry = 
a ia aa,’ OF dAy; Y dx; 
( dI dAy OF Ge) (48 — Ady) 
i=1 \Ol,’ dAé OAy; OAS Oxo; dAd 


(66) 
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In computing Aéd;;;, the moment 
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satis (= dl,’ OF dAy; OF =) 

i \Ol,dA5 Ay; dA6 AK, dAS 

(Ad — Ady_}) (67) 
So 
( 68) 
on ae ae " 
Aéby — Ad 

O< N N-1 <1 (69 


Aéy—; —Aéy_2 


As N increases, Ady approaches Aé, so this method jg 
convergent. 
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luids in Axially Symmetric Channels 


(Continued from page 156) 


and inlet ducts of turbomachinery. For instance, it is 
possible to determine the flow conditions in channels 
that are subdivided into several flow paths by means of 
surfaces of revolution with respect to the axis of sym- 
metry. . With proper modifications, it is possible to 
use this calculating principle for the determination of 
three-dimensional flow patterns in rotating wheels, 
either of compressors or turbines, where energy is added 


to or taken away from the fluid. 
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The Pressure on a Slender Body of 
Nonuniform Cross-Sectional Shape 
in Axial Supersonic Flow 


ERNEST W. GRAHAM# 
Douglas Aircraft Company, Inc. 


ABSTRACT 


The linearized theory for slender bodies is used to determine 
the pressure on an object that does not have rotational symmetry. 
The body extends infinitely far upstream as a circular cylinder, 
and in the region of interest (the distorted element) the cross 
section is gradually altered, without change of area, to a two 
lobed shape at the downstream end of the element. The series 
of shapes passed through in the transition from circular to two- 
lobed is fixed. However, the rate at which this transition occurs 
as a function of displacement in the downstream direction) may 
be varied. 

The formula for pressure coefficient shows the dependency on 
rate of change of shape and also indicates that the pressure co- 
efficient is proportional to the square of the thickness ratio of the 
distorted element. The thickness ratio of the distorted element 
is defined as the diameter of the initial circular cylinder divided 
by the length of the distorted element. 

Since the expressions used for the perturbation velocities are 
independent of Mach Number, the analysis applies to transonic 
and subsonic flows, as well as to the supersonic case. From this 
fact it can be shown that, within its limits of accuracy, the linear- 
ized theory for slender bodies predicts zero wave drag for dis- 
torted elements ending in cylinders. From reference 1 (and 
probably from reference 2) more general conclusions may be ob 
tained to the effect that altering only the cross-sectional shape of 
a body of revolution (and not the area) does not change the 
drag if the body ends in a cylinder or a point. From refer- 
ence 1 it can also be concluded that the dependency of pure shape 
distortion pressures on thickness ratio squared is general 

The pressure coefficients are calculated for a particular ex- 


ample and appear to be of reasonable magnitude. 


NOTATION 


x = displacement along the x-axis (in the flow direc- 
tion) 

l = length of the distorted element (in the flow direc- 
tion) 

g = (x/l) 

y = displacement along the y-axis 

g = displacement along the z-axis 

r = radius from the x-axis, = YY y? + 2? 

ro = radius of the initial circular cylinder 

R = (r/ro) 

6 angle of rotation about the x-axis 

U, = velocity at infinity (free stream) 

M. = Mach Number at infinity (free stream) 

\ = VM,?-1 

® = perturbation velocity potential 

Vir) = radial velocity (from x-axis) 

Veo) = tangential velocity (about the x-axis) 

Vix) = axial velocity (along the x-axis) 
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Vey) = axial velocity (along the y-axis) 

V(z) = axial velocity (along the z-axis) 

f(x) = a function representing the distribution of quadri- 
pole strength along the x-axis 

iL = (4/U.) fy*f(x) dx 

i* = (L/ro*) 

F(x, y, ) = a function determining body shape 

G = an arbitrary function arising in the determination 
of the body shape 

AW = perturbation of the absolute magnitude of the 
resultant velocity 

Ap = increment of static pressure from free-stream 
value 

q = dynamic pressure in free stream 

Cp = pressure coefficient = (Ap/q) 


Primes denote ordinary derivatives. Subscripts, unless other- 


wise indicated, denote partial derivatives. 


INTRODUCTION 


5 PURPOSE OF THIS STUDY is to determine, for 
slender bodies in axial flow, the order of magni- 
tude of the pressure and drag arising from lack of rota- 
tional symmetry. Two simple classes of bodies are 
suggested for investigation. 
that have the same cross-sectional areas at every sta- 
tion but not the same cross-sectional shape.* The 
second class consists of bodies having similar cross- 
sectional shapes at every station but varying cross- 
The first class is the subject of the 


The first consists of bodies 


sectional areas.‘ 
present study. 
For convenience, the type of distortion investigated 
is that produced by a simple quadripole distribution 
along the flow axis (a quadripole being produced from 
two dipoles as a dipole is produced from a source and 
It is considered that a circular cylinder extends 
The distribution 


sink). 
to infinity in the upstream direction. 
of quadripoles placed on the flow axis behind this cyl 
inder then gradually distorts the cylinder, producing, at 


some distance downstream, a two-lobed cross sec- 
tion. The cross-sectional area is not altered at any sta- 
tion. 


THE PERTURBATION VELOCITIES DUE TO A QUADRIPOLE 
DISTRIBUTION ON THE X-AXIS 
The linearized equation for perturbation velocity 
potential in Cartesian coordinates is 


[%,, + ®.,] — *@,, = 0 (] ) 
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Fic. 1. Coordinate system. 


or in cylindrical coordinates 


c ?, Pop ° 
&,, + + — | — r*,, = 0 (2) 
r 


r2 
where the expression in brackets is the Laplacian taken 
in a plane perpendicular to the flow direction and \ = 
V M,?—1. 

It has been indicated in reference 1 that slender body 
theory may be obtained by using solutions of the La- 
place equation in planes perpendicular to the flow. 
Certain restrictions are imposed on these solutions which 
correspond physically to a requirement that the body 
be smooth as well as slender. 

The expression 


® = —f(x) cos 26/r? (3) 


is a solution of the Laplace equation corresponding to a 
quadripole distribution on the x-axis, and a sufficient 














FLOW DIRECTION PERPENDICULAR 
TO THE PAPER 


Fic. 2. Cross sections of the distorted element. 


SCIENCES—MARCH, 1950 
restriction is that f(0) = f’(0) = f"(0) = 0 and tha 
x°f"(x) is not too large compared to f(x). [The mean. 
ing of these restrictions is made clear by substituting 
the expression for ® from Eq. (3) into Eq. (2) and noting 
that A°?,, in Eq. (2) should be small by comparison with 
other individual terms in the equation. | 
The corresponding perturbation velocities are 


®, =2f(x) cos 26/r° = V\,,) = radial velocity 

®,/7 = 2f(x) sin 20/r? = Vi = tangential velocity (4 
®, = —f'(x) cos 26/r?= V(,, —U,, =axial velocity 
DETERMINATION OF BoDy SHAPE FROM Bovunpary 


CONDITIONS 


In general, if F(x, y, 3) = O is the equation deter 
mining the body shape and V,, V(,), V(), ete., are the 
velocities in the indicated directions, then, since the 
velocity component perpendicular to the surface of the 
body must be zero, the following equation holds: 


OF _ 


OF V OF 
"Oz 


Vex) (y) 0 5 
Ox "Oy 

For the present problem, it is more convenient to use 
cylindrical coordinates and define the body shape by 


x= x(r, 0), giving 


: _. ‘Ox Vg Ox 
Vien — Vn =< 6 
Or r OO 
Letting I... = U., since the velocity perturbations 


are assumed small, and introducing the values of V;, 
and I’\) from Eqs. (4) gives 
2f(x) cos 200x 


r3 or r3 


2f(x) sin 200x : (7 
O00 
This equation may be linearized by the substitution 

L= (4/U.) J,* f (x) dx 
to give 


(2 cos 26/r°)L, + (2 sin 26/r4)L, —4 =0 (8 


The solution is obtained by solving the following equa- 


tions: 
r® dr/2 cos 20 = r4 d0/2 sin 26 = dL/4 (9 
The complete solution is 
L + (r* cos 20/sin? 20) = G(sin 26/r?) (10) 


where G is an arbitrary function. Assuming that the 
distortion is started from a circular cylinder of radius 
ry, then G can be determined from the condition that 
r = r when L = 0. 


(L/ro*) and R = 
(10) and solving for R, 


R -|) + 2L* cos 20+ V1 + 4L* cos 20+ = 


9 


It is convenient to let L* = 


(r/ro). Then, substituting into Eq. 


(11 
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PRESSURE ON A SLENDER 
This equation describes a particular cross-sectional 
shape for each value of L*, as shown in Fig. 2. If L* 
is chosen to be some function of x, a particular cross- 
sectional shape is produced at each station along the 
flow axis. 


DETERMINATION OF THE PRESSURE COEFFICIENT 


If AW (the perturbation of the absolute magnitude 
of the resultant velocity) is small by comparison with 


the free stream velocity, U’.., then 


(12) 


C, = —(2AW/U.) 


Retaining the radial and tangential velocities affect 
ing AW gives, with sufficient accuracy, 


ey ey] 
U. U. rU. , 


Substituting into Eq. (13) the perturbation velocities 
corresponding to a quadripole distribution along the 


C, = 


x-axis [see Eqs. (4)], 


C. = 2f’(x) cos 20 (A) 
p> r2U rU. 


(14) 


For convenience, let ] be the length of the distorted 
eement and let € = (x/l). Since f(x) = (U,/4) 
(dL/dx) with L = L*ryi and since R = 7/71, the follow 
ing equation is obtained from Eq. (14): 


(“) i cos 26 (— -) l | (3 
_ 5) 
l dé? 2R? dé J 4R°® 


If (27)/l) is defined as the thickness ratio of the dis- 


C, = 


torted element, then the pressure coefficient is propor- 
tional to the square of the thickness ratio. Although 
Eq. (15) applies only to the specific type of cross-sec- 
tional shape distortion considered here, the dependency 
of pressures due to shape distortion on the square of the 
thickness ratio is actually general, as can be demon- 
strated from the equations in reference 1. 


From Eq. (15) it is possible to show in some specific 
cases that the distorted elements have zero wave drag. 
However, from Eqs. (4) it is seen that the perturbation 
velocities are independent of Mach Number, and the 
analysis therefore applies to transonic and subsonic 
flows, as well as to the supersonic case. It is known that 
a body ending in a cylinder has no pressure drag in an 
incompressible flow. Then, since the present method 
predicts the same pressure drag for supersonic as for 
incompressible flow, it cannot accurately predict any 
pressure drag supersonically for bodies ending in a 
cylinder. Reference 1 shows that, for a body of revolu- 
tion ending in a cylinder or a point, the drag is not 
altered by changing only the cross-sectional shape. 
This is true within the accuracy of the 
linearized theory for slender bodies. An equation 
for drag obtained by Ward? suggests a similar re- 
sult. 


limits of 














BODY IN SUPERSONIC FLOW 175 
CIRCULAR 
—* CuLINDERm)2— DISTORTED ELEMENT——~— 
mi wes — 
| 2 4 
x r 
e o® 






~~ 

~ 

“ 
"Stlatinenaet 





s+! 








FLOW DIRECTION 
Shape of element (example with radius magnified 


Fic. 3. 


EXAMPLE 


As an example, consider a body that starts as a circu- 
lar cylinder and undergoes a gradual transition to a two- 
lobed cylinder, as defined by Eq. (16): 


L* = (1/4)}1 — cos [wé2(2E — 3)]} (16) 
The distorted element extends from ¢ = 0 to € = 1, 
since L* = OQ até = Oand L* = 1/2at&é=1. Itcan 


be checked by differentiation that 
L*'(0) = L**(0) = L*’’"(0) = 0 


which was specified as a requirement for smoothness. 


(Continued on page 192) 


















































10 
5 
= ° RS 
fe /A) @:2 ry, 
015 arm 
Dp O= o? 
-§ 
1.5 
=0° 
g 245° 
1.0 
9 
R= <9 
6 4 
0.5 
oO 
0.25 0.50 0.75 1.00 
§=(*#) 
Fic. 4. Ordinates and pressures (example). 








Reply to Discussions of the Wright Brothers 
Lecture’ 


A. E. RUSSELLt 
The Bristol Aeroplane Company, Ltd. 


A NUMBER OF THE POINTS raised in the discussion 
are concerned with the inherent restrictions on 
the operation of turbine-powered aircraft. As Mr. 
Kelly comments, if for any reason it becomes necessary 
to depart from the flight plan, there is, by comparison 
with piston-engined aeroplanes, much reduced flexibil- 
ity of choice of operating speed and altitude. Mr. 
Hazen points out that some alleviation of this handicap 
is possible on a multiengined aeroplane by shutting 
down engines and, that with a coupled installation, 
this can be done without stopping propellers; so the 
need for an apology to the passengers would not arise. 

For an aeroplane with eight engines, for each engine 
shut down the best altitude falls by about 4,000 ft. at 
the expense of about 25 m.p.h. in cruising speed and 
an increase of about 4 per cent of fuel consumed. In 
any case, reserves must be carried to provide for the 
possibility of inadvertent engine failure. Whereas at 
first sight the procedure of selecting the best combina- 
tion may appear to place an added burden on the flight 
crew, if ambient conditions are known, there is no 
ambiguity as to the right course of action. 

The argument Mr. Kelly used in his comparison of 
the net performance of aeroplanes powered with piston, 
turboprop, and turbojet engines is appropriate and 
can be extended to the case of greater ranges. It 
emphasises the danger of too general a conclusion and 
the increasing need for tailoring an aeroplane for specific 
routes. For example, if intermediate stops are com- 
pulsory for refuelling, then, in addition to the effect 
on block speed, time is lost on the ground and the cruis- 
ing speed at altitude is a quite unreliable measure for 
journey time. 

Also, the daily and seasonal variation of wind speeds, 
known to increase in velocity with altitude, may in 
some cases disproportionately penalise turbine aircraft 
with fuel reserves. The less flexible aeroplane suffers 
by its reduced ability to overcome adverse atmospheric 
conditions. 

In making a case for the use of turboprops for short 
range medium-altitude operation, Mr. Fink calls for a 
low-compression engine. It appears likely that, in 
selecting the best engine characteristics, the value of the 
compression ratio should bear some direct relationship 





* “Some Factors Affecting Large Transport Aeroplanes with 
Turboprop Engines,”’ by A. E. Russell, JouRNAL, Vol. 17, No. 2, 
pp. 67-106, February, 1950. 

+ Chief Designer, Aircraft Division. 
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to the take-off power. 
units in pairs, a relatively small range of basic engine 


As it is convenient to coupk 


types will cover a large range of aeroplane requirements 

Mr. Hazen suggests that aircraft development may 
not be able to keep pace with the rapid rate of develop. 
ment of turbine engines. Firstly, it will never be 
difficult to persuade an aircraft man to accept a lower 
fuel consumption. If improved consumption is also 
accompanied by greater power, then, if there is a tem- 
porary limit on gross weight, the mean cruising altitude 
is raised. A general idea of the changed conditions 
can be had by inspection of the curves given for the 
engine inoperative case, but here, of course, the changes 
are of opposite sign. Higher power, besides increasing 
the altitude, raises the cruising speed and slightly 
improves fuel consumption. If toward the end of the 
flight the altitude must be limited for structural reasons 
the cost in fuel can be estimated from the cruising 
chart. 

It is unlikely that power-plant development would so 
outstrip aircraft development that even with a multi- 
engined arrangement it would be feasible or desirable 
to eliminate a power section. 

On another point, I am sorry if I have misled anyone 
into supposing that the effect of flexibility of large 
aeroplanes at high speed will have any major effect 
on passenger comfort. The effect is restricted to the 
relative reduction of the amplitude of wing flexure due 
to aerodynamic damping which has a much larger 
influence on wing bending moment than on vertical 
acceleration at the centre of gravity. 

Mr. Silverstein suggests that it is difficult to justify 
the application of turboprops for long range with the 
modest engine characteristics quoted in the lecture. 
However, starting with power-plant performance im- 
mediately available enables a general picture to be 
drawn, in which it is possible to place into some pet- 
spective the various related but individual problems. 
With so much acknowledged scope ahead for the devel- 
opment of turboprops, we can gain confidence that 
this form of power plant will produce a good balance 
between high performance, comfort, and operating 
economy over a wide range of stage distances. While 
a strong case can be made for other forms of power 
plants for the same and/or other operations, it is dis- 
turbing to reflect that, with a finite capacity for devel- 
opment effort, a dispersal of efforts retards development 
in each direction. 
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REPLY TO D 


The problem of gusts is in a different category, for 
here our concern is confined to refinement in marginal 
At this stage we must question the 


circumstances. 
I agree 


validity of both premise and form of analysis. 
with Mr. Silverstein that detectors on the fuselage 
will recognise rogue, as well as formalised, gusts. In 
the end, as with other aspects of safety, we must fall 
back on probabilities, but as yet we lack authentic 
statistics. 

The use of the gust alleviator is not so critical as 
Professor Hoff suggests. Under the standard maxi- 
mum gust for the aeroplane in question, the factor of 
safety exceeds unity at all speeds. The maximum 
potential alleviation greatly exceeds the needs of nom- 
inal strength requirements, so that a result that shows 
about 30 per cent effectiveness is adequate. The task 
of simplifying the signalling system is not discouraging; 
oue might add that here there is plenty of scope. 

On the question of fuselage shell stability, we are 
indebted to the work of Professor Hoff. If the frames 
appear to him to be too close, their position was dic 
tated by window and door locations. The ultimate 
strength was finally established by a test to destruction 
on a half-scale model. 

I believe that, since the example mentioned by 
Professor Hoff, much has been done to direct the trans- 
between our two 

Hitherto, it was 
The most effective 


mission of scientific information 
countries into accessible channels. 
frequently filed in obscure archives. 
contact is always personal, and the now established 
Joint Conferences should do much to reinforce the 
official procedure. 

I find myself in accord with Mr. Wells’ remarks, for, 
when it is necessary to introduce time as an additional 
parameter into structural analysis, we are ill prepared 
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to assess the appropriate loading history. We are 
merely aware of the phenomenon but without means of 
assessing an empiricism. It is not essential for a set 
of design calculations to assume a close approach to 
precision, but it is important that departures from 
convention, in aircraft size or shape and in operating 
circumstances, shall have proper consideration on a 
logical basis. I share Mr. Wells’ concern for the need 
to establish design data for gust loading for both sym- 
metric and asymmetric application. 

It is certainly fortunate that sweptback wings have 
built-in gust alleviation. It is also unusual for any 
arrangement dictated essentially by aerodynamic causes 
to have beneficial structural effects. These are out 
lined concisely in the discussion. The magnitude of 
such alleviation permits more than local refinement, 
for it plays a part in the initial selection of the design 
variables. 

I am very encouraged to hear that the testing of 
scale flutter models is an effective design aid. Clearly, 
any form of calculation is limited by the validity of 
basic data; hence, if it can be argued that an experi- 
mental method reasonably simulates actual conditions 
or even if it merely indicates qualitative effects of 
changes, then the advantage of visual illustration and 
speed of assessment is obvious. The examples given 
by Mr. Wells not only justify scale flutter models but 
place them in the top priority of design methods. 

I wish to thank all those who have contributed to 
this discussion, for their personal comments, and for 
the way in which they have amplified the subject matter. 
If I have succeeded in promoting a basis for further 
understanding, then I am proud of a small share in 
the unceasing flow of information, eastwards and west- 
wards, across the Atlantic. 
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Station Functions in Flutter Analysis 


J. H. Greidanus and A. |. van de Vooren 
National Aeronautical Research Institute, Amsterdam 
October 15, 1949 


- THE FIRST PART of the paper ‘Station Functions and Air 
Density Variations in Flutter Analysis’? (JOURNAL OF THE 
AERONAUTICAL SCIENCES, Vol. 16, pp. 345-353, June, 1949), 
Manfred Rauscher has presented a method of calculating fre- 
quencies and modes of continuous systems of which the essential 
features are: 

(i) The determination of an interpolation function f(y) to the 
deflections 2;(4 = 1, 2,..., ) at a series of reference stations y;, 
representable by 

n 
f(y) = yi sifi(y), OS 951 (1) 
i=1 
where the fi(y) are station functions—i.e., functions which vanish 
at all reference stations except one, this station being different 
for every function of the set. 

(ii) The identification of assumed and calculated deflections 
at the reference stations, the calculated deflections being obtained 
from inertia loads based on the assumption (1) by means of in- 
fluence functions. 

It is suggested, though not stated explicitly, that the reference 
stations should divide the interval in equal parts. It is recom- 
mended to identify the station functions with continuous func 
tions, having continuous derivatives up to the second order (in 
the case of bending), satisfying the geometrical and one of the 
dynamic boundary conditions, and representable within each 
subinterval by polynomials in y of lowest possible degree. 

Now the polynomial approximation of modes is a classical pro 
cedure of vibration analysis."? Using one polynomial in the 
full interval, imposing the full set of geometrical and dynamic 
boundary conditions as usual, and restricting the degree (again 
for bending vibrations) to the fourth, fifth, etc., for subsequent 
exemplars from the required set, the following well-known three 
term polynomials result :! 


1/6(n + 2) (n + 3)y"t! — 1/3n(n + 3)y"*t? + 
1/6n(n + 1)y"*3; n= 1,2, etc. (2) 


Comparing this method with Rauschers’, it is not clear what 
kind of advantages should result from the neglect of one dynamic 
boundary condition and the introduction of discontinuities, at 
the reference stations, in the third-order derivative. The classi- 
cal polynomials do not have the station-function characteristic, 
but could easily be made to do so by suitable linear recombination. 
There is, in general, no advantage in this reduction, but Rauscher 
states that gains in convenience result, connected with the rep- 
resentation of unspecified aerodynamic finite. span corrections in 


flutter calculations. Now, since these corrections, to be applie 
to the aerodynamic loading, will doubtlessly be linear, and since 
they are stated to vanish with tip deflection, they are evident) 
proportional to tip deflection, the probably complex ratio de 
pending on y (giving the spanwise distribution of the correctig, 
and perhaps on reduced frequency and aspect ratio. If so, y 
labor saving of any importance seems imaginable. 


The proposition to improve station functions by application g 
the transformation 


1 
fin) = f K(n, y)m(y)fi© (y)dy 


where A(n, y) = influence function, m(y) = mass distribution 
agrees again with a well-known procedure of vibration analysis, 
applicable with equivalent success to all ‘‘zero-order”’ approxi. 
mations of modes. The transformation reduces the components 
of disturbing modes and leads, when iterated, to the fundamenta| 
mode. When not only the fundamental but also overtones mys 
be calculated, its application in combination with collocation or 
Rayleigh type analysis should be handled with some care, sing 
otherwise even desired modes may disappear. To prevent 
this, computations should be made with an increased number of 
figures. 


Concerning recommendation (ii), it is not clear why the colloca 
tion method has been preferred to the Rayleigh-Ritz-Galerkin 
method, which has the beneficial property of second-order errors 
in the frequencies, when the errors in the modes are considered to 
be of the first order, at least if the differential equation (s) of the 
system is (are) self-adjoint. The recommended method generally 
leads to first-order frequency errors. 


When for some reason or other, use of collocation method is 
preferred, it can be shown to yield the best results when the collo 
cation points coincide with the nodes of the first disturbing mode 
In the example of the uniform beam of which the lowest two ben¢- 
ing frequencies and modes have been calculated, the collocation 
points should be taken at the nodes of the third mode—i.., at 
y ~ 0.50 and y ~ 0.87. Though one should be very careful in 
checking the theoretical considerations developed above by one 
single example, it may be of interest to compare results for the 





calculation of the lowest two bending frequencies 1, v2 of the 
uniform beam. 
Vv; yp» 
I 3.52197 22 . 2789 
II 3.51600 22.1780 
ill 3.51721 22.1440 
IV 3.51604 22.7127 
V 3.51602 22 .0345 
I = station functions, collocation at y = 0.5 and y = l 
(Rvuscher); II = same station functions, collocation at y = 0 
and y = 087; III = same station functions, Rayleigh; IV = 
polynomials (2) form = 1 and m = 2, Rayleigh; V = exact 


All results are as expected, except that the continuous poly 
nomials appear to be less suitable for an accurate representation 
of the overtone. This feature should, however, be considered 4 
purely incidental and the same polynomials will doubtlessly be 
particularly efficient for certain beams of different stiffness of 
mass distribution. Finally, it may be remarked that, though a0 
analysis on concentrated-mass basis requires a far greater num- 
ber of reference sections, useless roots never are calculated whet 
applying a matrix iteration method instead of solving a chat 
acteristic equation. Performed in this way, a concentrated 
mass analysis remains a possible alternative to the Rayleigh 
analysis with simple polynomials. 


REFERENCES 


1 Duncan, W. J., Galerkin’s Method in Mechanics and Differential Eq 
tions, R. and M. 1798, 1937 
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READERS’ 


:Biezeno, C. B., and Grammel, R., Technische Dynamik, p. 747; 


|, Springer, Berlin, 1939 
" 4 Biezeno, C. B., and Grammel, R., loc. cit., p. 169 


Author's Reply 

In Part I of my paper on ‘Station Functions and Air Density 
Variations in Flutter Analysis,” to which the preceding note 
refers, a relatively large amount of space is given to the develop- 
ment of the interpolation idea from elementary beginnings. 
The example of a uniform cantilever beam is used to show how 
the idea works out under simple conditions and how the accuracy 
of the outcome depends on the care of the interpolation. It is 
then observed that the use of standard polynomials, while 
atisfactory with the continuous beam, will prove unsatisfactory 
when there are discontinuities both in the elastic and in the 
inertial properties of a system—hence, the suggestion that the 
interpolation functions be adapted to the peculiarities of the 
stem. This adaptation proposal, not the passing glance at 
the familiar polynomials, is the cardinal feature of the discussion. 

A careful reader of the original paper might have noted that 
the polynomials chosen for illustration were of third degree, 
instead of fourth or higher, because they were intended to be 
intermediate between the straight-line functions on the one 
hand (which did obvious and serious violence to the boundary 
conditions) and the functions based on triangular loadings on 
the other hand (which took the boundary conditions into me- 
ticulous account). There was no expectation, implied or ex- 
pressed, that any advantages, beyond a simplification of the 
interpolation functions, should result from the neglect of one of 
thedynamic boundary conditions in this case. As the condition 
violated was that of continuity of shear, the inaccuracies resulting 
from the violation are traceable to fictitious concentrations of 
lads at the points concerned. The quality of the results thus 
permits a judging of the possible satisfactoriness of an experi 
mental determination of interpolation functions by means of 
concentrated loadings on a structure. Such a determination 
has been used with success in a Station Function analysis of the 
vibrations of a sweptback plate.' 

With regard to the recommended procedure of adapting the 
Station Functions to the nature of a system, Greidanus and 
van de Vooren take the view that the scheme is simply the one 
underlying the well-known transformation formula they quote 
A more careful scrutiny of the proposed plan should make clear 
to them that this view is mistaken. The determination of the 
functions in the manner actually suggested in the paper involves 
ueither collocation nor any other type of explicit comparison 
between a first and a second set of deflections. The subsequent 
use of the functions in a collocation operation according to the 
formula quoted does not serve to supply or improve the Station 
Functions. Iterations of 
finement, with all their familiar troubles of mode suppression, 
tte, are deliberately avoided through the painstaking initial 


and other processes successive re- 


choice of the functions. 

The reduction of interpolation functions of any sort to the 
form of Station Functions was admitted to be a matter of pref- 
fence, not necessity. The Station form of the ex- 
pressions does, however, have the advantages of maximum 


Function 


Physical clarity and of allowing the most convenient tie-in be- 
tween the motions of the distributed masses and of masses con 
centrated at the reference stations. This tie-in makes it natural 
to choose the stations at points of mass concentration. It is 
emphatically not suggested, as Greidanus and van de Vooren 
lake it to be, that the reference stations should always have 
intervals them 
usually be chosen when there are no reasons for another pref 


uniform between though equal spacings will 
erence, 
The methods of Rayleigh-Ritz, Galerkin, and Least Squares 


were not overlooked by me. Twenty years ago,” I suggested a 
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simple Rayleigh-type procedure that does not seem too hopelessly 
out of date today. Concrete experience, not abstract specula- 
tion, has led to the belief that the Station Function method offers 
an especially favorable combination of convenience and accuracy 
for the unified handling of the still-air and moving-air phases of 
the problem. As the interpolation curves adjusted to a system’s 
properties already are close to the exact deflection curves, there 
would actually be little gain from an averaging out of the differ- 
ences by an integration process instead of the proposed method 
of simple collocation. The table of comparative results in the 
preceding note makes this clear even with the crude Station 
Functions that are taken for illustration. 

It is recognized, of course, that the determination of the sug- 
gested interpolation curves by integration of static load-deflection 
curves requires a substantial advance investment of labor. This 
investment is, however, smaller than that commonly made to 
determine interpolation curves in the form of natural mode 
lines; it is believed to result in a gain in economy, as well as a 
gain of transparency, through the unification of method it 
permits in the handling of the flutter problem as one whole 

Admittedly, the discrete-point method offers an alternative 
way of treating the still-air and the moving-air phases of the 
problem according to a single general plan. This alternative may 
remain appealing to someone so accustomed to it that he no 
longer notices its tediousness. But the fact is that the numerous 
useless roots, even though they need not be calculated, do impede 
the solution for the useful roots. Also, the matrix iteration 
process frequently is held up by the proximity of two roots, 
especially those corresponding to the second and third coupled 
modes of a typical wing, and the whole iteration plan bogs down 
discouragingly when the analysis advances from the still-air to 
the moving-air phase of the problem. Still, in spite of these 
shortcomings, there certainly is no denying of the method’s 
record of great utility. 

Corrections of the aerodynamic effects for the influences of 
finite span have been proposed by various authors, notably 
Reissner, Biot, Jones, Kiissner, and Cicala. My research group 
has been in closest touch with Reissner and is most familiar with 
his method of correction.* The finite span effects are being 
considered with increasing seriousness in advanced American 
practice. Because their incorporation into the analysis involves 
considerable extra labor, a careful planning of the procedure is 
of some importance. The desirability of having all but one of 
the deflection functions go to zero at the tip, which Greidanus 
and van de Vooren question, comes from the fact that the Fourier 
series for all but one of the component circulation or load dis- 
tributions will then be rapidly convergent. It is true, of course, 
that the aerodynamic loading must come to the same total, na 
matter how it is built up of component loadings. The benefit 
of the plan advocated in the original paper is that it does simplify 
the computation work. 

As far as the overall worth of the Station Function method is 
concerned, a fair judgment can naturally not be formed on the 


For this reason, the 


scant basis of an outline description. 
original paper made clear reference to the fuller reports submitted 
to the Bureau of Aeronautics, which also contain a fuller bibliog- 
raphy. A consultation of the reports in question and direct 
communication with me and my associates at M.I.T. are sug 
gested as the best ways of clearing up points that may still remain 
obscure to anyone desirous of getting all the known facts about 


the method 
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A Note on the Velocity of Sound 


M. V. Morkovin 
University of Michigan, Ann Arbor, Mich. 
October 17, 1949 


EFERENCE | called attention to a certain lack of rigor found 

in most of the elementary derivations of the formula for 
the speed of sound. The main point of the criticism was that 
the ratio dp/dp was indeterminate in the sound wave and should 
not be treated as a derivative. While agreeing about the lack of 
rigor in the customary derivations, the present note concludes 
that the key to the confusion is not the alleged indeterminateness 
of dp/dp but the lack of clarity in specifying the variables in- 
volved, especially the independent variables. In his counter- 
example, the author of reference 1 unfortunately forgets which 
quantities are variable and which are constant under the original 
assumptions of his problem. 

(1) The assumptions of the problem are those of steady, 
continuous, adiabatic flow, without viscosity, in a duct of constant 
area. Then and only then Eggs. (1), (2), and (3) of reference 1 
namely, p du + udp = 0,pudu+dp =0,andudu + C,dT = 
O—are valid simultaneously. These assumptions imply a flow 
at constant entropy, s, as the elimination of the term C, dT 
[which equals T ds + (dp/p)] between Eqs. (2) and (3) demon- 
strates. As the author of reference 1 remarks, these equations 
admit only the “‘trivial’” solution « = constant, p = constant, 
p = constant; in an isentropic flow through a one-dimensional 
duct, there is simply no agency that could cause a change in 
these quantities. 

The counterexample of reference 1 considers the particular 
solution of these equations, u2 = yp/p. This is one of the 
family of the constant solutions above, for which the constant 
values of the quantities u, p, and p are numerically related. 
Under the assumptions of the duct problem du = 0, dp = 0, 
and dp = 0, Eqs. (1), (2), and (3) are satisfied without any contra- 
diction. 

In order to allow a continuous change of u, p, and p, in steady 
flow through a constant area duct, as reference 1 assumes without 
explicit statement, one of the other assumptions must yield. 
The usage of Eq. (2) in the counterexample implies that the 
that heat 

yp/p as an identity for 


flow is still assumed nonviscous but may be added 


steadily. This and the usage of u? 
all flows in the given duct of necessity leads to nonisentropic 
processes characterized by Eq. (7) of reference 1. There is no 
coutradiction, merely a misunderstanding. 

(2) When more complicated flows are considered, such as 
viscous and diabatic flows, a more careful specification of the 
thermodynamic process corresponding to the phenomenon of 
propagation of sound is in order. The properties of a homog 
enous medium that is in thermodynamic equilibrium depend 
upon two independent variables, say p and s. In particular, 
Pp = f(p, s). Following reference 2, a complete and rigorous 


definition of the speed of sound is obtained through the identity 


°*) 
a? fo(p, 5) (4) 


It is immediately recognized that the use of the ratio dp/dp 
from Eqs. (1), (2), and (3) is covered by definition (4). 

(3) In order to dispose of the question of indeterminateness 
of the ratio dp/dp, it is merely necessary to use exactly the same 
reasoning as reference 1 did in its derivation of the limit 1,2 —- 
yp:/p1 as a shock wave approaches a sound wave. A somewhat 
tedious manipulation of the three conservation laws valid across 
a shock wave? leads to 

Pi-pi _ pi + [lp — p)/2) (5) 
2— ~\ pi + [(o2 — p1)/2] 


As the strength of the shock wave vanishes, the ratio in question 


AERONAUTICAL 


SCIENCES—-MARCH, 1950 


approaches not only a determinate value but also the “correct 
one, ypi/p1-. 
concept of a sound wave as a limit of a shock wave checks with 
definition (4). Incidentally, Eq. (5) is unsteady 
flows, and our deductions are general. 

(4) Finally, it is proper to remark that Lamb’s “rigoroy; 
derivation of the equation for the speed of sound, recommende 


Since, furthermore, (s: — s,) approaches ZeTO, the 


valid for 


by reference 1, is subject to exactly the same criticism as th 
elementary derivations originally criticized. 
determinateness of dp/dp and a specific thermodynamic proves 
with p as independent variable when he replaces dp/dx by 
(dp/dp)(Op/0x). 
sure and density changes identifies the process as isentropic 


Lamb‘ assumes the 


The further assumption of infinitesimal pres 
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Further Remarks on the Velocity of Sound 


Morris Morduchow 
Department of Aeronautical Engineering and Applied Mechanic 
Polytechnic Institute of Brooklyn, Brooklyn, N.Y 


November, 1949 


I SHOULD LIKE, if possible, to eliminate any confusion that 

may arise from a comparison of references 1 and 2. 

In reference 1, item (1), it is at first agreed that Eqs. (1), (2 
and (3) of reference 2 imply u, p, p, and T to be constants (wit 
the space variable x). This, however, appears to be forgotter 
in the second part of item (1), reference 1, where it is attempte 
to explain the contradiction pointed out by Eqs. (4)—(7) inréd 
erence 2 with the statement that, if the equation 


u? = y(p/p) 
is treated as an identity, then this merely replaces Eq. (3) 0 
The fact 
that is overlooked in reference 1 is that the set of Eqs. (1) ant 
(2) of reference 2, and Eq. (a) will also imply only the “trivial 


reference 1 and implies a ‘‘nonisentropic process.’ 


type of solution: Consequently, these equa 
tions cannot imply any “‘process” at all in which p may & 


It is implied in re 


p, ep, u constant. 


cousidered as varying as a function of p. 
erence 1 that Eq. (7) of reference 2 represents such a “‘process 
However, a different and contradictory “‘process’’ could hay 
been mathematically arrived at by substituting wu du 

(1/2)yd(p/p) into Eq. (1), instead of Eq. (2), of reference - 
Then one would obtain p = k/p, contradicting Eq. (7). Mor 
over, contrary to the belief in reference 1, even an ‘‘isentropt 
could be mathematically arrived at by substitutimg 
(1) ane 


process”’ 
Eq. (a) into Eq. (4) 
(2)] of reference 2. As implied in reference 2, such apparel! 


[which was obtained from Eqs 


contradictions are bound to occur when equations like E® 
(1), (2), and (3) of reference 2—or Eqs. (1), (2), and (a)— 
incorrectly treated as if the quantities p and wu could be co 
sidered as dependent variables with p as independent variable 
instead of remembering that all of the quantities u, p, 7, - 
p are really constants in these equations, with the space variable 
x as the independent variable. (If there does exist any confusio 
with respect to dependent and independent variables, it is thus 


in reference 1.) 
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In reference 1, it is agreed that ‘‘elementary derivations” often 
lack rigor, but the author is fairly vague with respect to just 
where the rigor is lacking. Perhaps the entire discussion can be 
clarified and all apparent inconsistencies avoided if it is realized 
that a set of equations like Eqs. (1), (2), and (3) of reference 2 
can be readily transformed into a set of algebraic equations in 
which all quantities are simply constants. This can be done 
eg. reference 2) by establishing the shock-wave equations 
and then considering a sound wave as the limit of a weak shock 
wave. Reference 1 [item (3)] is in agreement with reference 2 
in this regard. However, since all quantities are identically 
constants, it is doubtful that differential expressions such as 
qt = dp/dp can be properly derived and interpreted by such 
means. 

With respect to item (4) of reference 1, it may be remarked 
that the only assumption in the derivation of the sound velocity 
inreference 3 is that p is a function of p and of ponly. Although 
such an assumption, as shown in reference 2 and in the present 
note, leads to mathematical inconsistencies in the steady one- 
dimensional differential equations of nonviscous flow without 
heat addition, it does not lead to such inconsistencies in the 
unsteady-flow equations, since the latter equations admit of 
general solutions other than the “‘trivial’”’ type. 
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On the Compressibility Correction for 
Subsonic Unsteady Flow 


John W. Miles 
Department of Engineering, University of California at Los Angeles 
September 26, 1949 


I A PREVIOUS ITEM in this column,' we set forth a transforma- 
tion that reduced the problem of low-frequency harmonic 
motion of a supersonic airfoil to an equivalent problem in steady 
flow. A general study of low-frequency unsteady-flow effects 
at all speeds was recently (Summer, 1949) undertaken at the 
Aerophysics Laboratory of North American Aviation, Inc., and 
asa part of this study, a similar transformation has been applied 
to reduce low-frequency unsteady-flow problems in sub 
sonic compressible flow to equivalent problems in incompressible 


flow 
The linearized boundary value problem for the oscillating wing 
may be pesed as 


a 2 
Vv? — ur ( + it) o(x,y,z) = 0 (1) 
Ox 


$:(x,y,0+) = w(x,y), (x,y) in S; (2) 


Ox,9,0 +) = +(x1,y,0+) exp[—ik(x — x)], (x,y) in Sp (3) 
o(x,y,0) = 0, (x,y) in S (4) 

k = (wb/U) (5 

where a flow of velocity U and Mach Number M is directed 
along the x axis over an airfoil for which the projection on z = 0 


is S, and which has the prescribed normal velocity w(x,y) exp- 
The regions S, and S; represent the wake [including the 
xi(y)] aft of the wing (S,) and the remainder 


tt), 


trailing edge, x = 
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of the plane s = 0, respectively; & is the reduced frequency; and 
b is the unit of length. Eq. (3) states that the potential is dis- 
continuous across the vortex sheet constituting the wake, while 


the pressure is continuous there. If we introduce the transforma- 


tion 

x’ =x, y’ = By, 2’ = Bs, b’'=B-'b (6) 
o'(x’,y’,2’) = exp(—1dAx)o(x,y,2) (7) 
w’(x’,y’) = B-! exp(—?Ax)w(x,y) (8) 
k’ = (k/B?) (9) 
B= (1 — M?)' (10) 

we find that Eq. (1) transforms to 
[V2 + (k’M)?]o'(x’,y’,2') = 0 (11) 


while the boundary conditions (2), (3), and (4) remain invariant, 
except for the scale transformation of S;, S:, and S;. Hence, if 
we neglect terms of order (k’ 17)? (and possibly, order k? In k), our 
problem is reduced to one in incompressible flow. The trans- 


formation law for the potential follows as 


att He ae sae k ut = B-! [1 + 7A(x/b)] X 

lb’ sb’ 6b’ US f ; a 

xy 3 _ w aes ee 
\ b’ Bsb’ Bb’ Ben U’ (4/6%), +e Ra 


(12) 

s and 6 are aspect and fineness ratios, respectively. In general, 

the dependence of ¢ on (w/U) is linear and homogeneous, the no- 

tation being operational. The transformation law for the pres- 
sure on the wing is found to be 

ix ¥ w fx 

p < . k, > 1 p d 

Lb’ asd’ I \ 


ae ee 
(k/B*), O> 
}’ Bsb’ Uv’ ; { 


x 9 Ww 
(x/o)p ‘ A oF ooh 


Bsb lL 
ix ¥ ( “) ' . 
. ° 2 Os 
P to’ 8sb Nb OU f 


Our extension of the Prandtl-Glauert rule, as represented by 


i (RAL?/B%) 


+ 0 [(kAl/B)2] (13) 


Eq. (13), may now be stated as follows: To find the pressure dis- 
tribution on an airfoil of aspect ratio s executing harmonic mo- 
tion at reduced frequency & in a subsonic compressible flow, we 
calculate the pressure distribution on a second airfoil of aspect 
ratio Bs and reduced frequency k/8? executing the same motion, 
add to this the distribution on the same airfoil multiplied by 
iA(x/b), subtract the distribution for the second airfoil with a 
new downwash 7A(x/b)w, neglect unsteady flow effects in the last 
two distributions, and divide the overall result by 8. In carrying 
out these operations, it should be specifically remarked that any 
frequency dependence exhibited by w is not to be modified 
throughout the transformation, except insofar as parts of it may 
be neglected in conformance with the previous approximations. 

It should be understood that these above results are not ap 
plicable to two-dimensional problems, since the solution for 
would involve integration over the infinite span, and this integra- 
tion generally would not be commutative with the assumed ex 
pansion in powers of the frequency. The solution for low- 
frequency compressible flow has been given earlier,? and the 
compressibility correction obtained therein is rather more com- 
plex than that outlined above. Moreover, it is of particular 
interest that the logarithmic frequency term entering the two- 
dimensional solution is anomalous, and an examination of the 
low-frequency limit of the three-dimensional results? shows this 
logarithmic term to be important only for very large aspect 
ratios (say, s > 10). 

Finally, we remark that the foregoing results are valid only for 
planar (wing) problems because of the assumed form of the 


boundary conditions. 
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Upper and Lower Bounds on Variational 
Integral in Compressible Flow 


G. V. R. Rao 
New York University, New York, N.Y. 
November 1, 1949 


ECAUSE OF THE INHERENT DIFFICULTIES in obtaining solu- 
tions for compressible flow equations, we have no way of 
estimating the error in the solution. By variational methods,* 
an approximate solution can be obtained. 
an integral, the first variation of which leads to a linear differ- 
ential equation for the perturbation potential. When the do- 
main extends to infinity, the integral has to be properly corrected 
by the addition of a boundary integral as indicated by Wang.’ 
The correct form of the variational integral would then be as 
A value of y = 2 is used for simplicity in represen- 


OF d¢i Op» \? 
r ‘ ) - . [( : ) r 
Oy Oy Ox 


06 oe)" | 


Braun! formulated 


given below. 
ation. 


t 
f (2 Odi 
I, = 4 Qh, 
cof fame 
Od» \? Odi Ode 
(s)]ee(e 
oO”, Ox OX 


( dx dy (1) 


oy Oy 
where 
¢: = velocity potential for incompressible flow 
¢2 = perturbation potential 
F = qm? — (0¢i/0x) — (O¢i/Oy)? 
gm = Maximum velocity that could be obtained by expansion 
into vacuum 
The approximation involved in using y = 2 was investigated, 


by Wang.* Because the present problem is linear in its nature, 
a value of y = 1.5 or 1'/; could be used without involving elabo- 
rate computations. 

By the application of Friedrichs’ method? to the above in- 
tegral, an adjoint minimum problem is formulated. The inte 
grand of the above variational integral is denoted by G, and we 


set 


OG Ou. = a; OG OV» = 5 (2) 


where #2 = O¢2/Ox and v2 = O¢2/dy. 
It can be easily verified that the first variation of the integral 
leads to the differential equation 


(O0a/Ox) + (0b/dy) — f = 0 (3) 
where 
fb (°F do; OF 22) 
" “Now ox oy Oy 
The natural boundary condition 0¢./dn = 0 on the body can 
be expressed as 
a(Ox/On) + b(dy/dn) = O (4) 


The relations (2) can be solved explicitly for #2 and 2» 
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The conjugate functional can be written as 


ICAL MARCH, 


“i = 


and 


Jt = (1/2)(usa + v2b)dx dy f 


with #2 and v2 expressed explicitly in terms of a and b. We wi 
choose 
a = 2F (0¢;/0x) + py ) : 
b = 2F (0¢i/dy) — yz J 
By this artifice the integral (( 


satisfying Eq. (3), a priori. 


can be expressed in terms of derivatives of y. For the functio; 
¥Y, we can choose any function of order less than 1/r, and satis 
fying the boundary condition dy /0s = O on the body. 

Thus we have a new variational integral (6) which, together 
with the integral (1), gives us bounds on the value of the integral 
From a solution obtained by minimizing integral (6), we can 
with the aid of Eqs. (5) and (7), gather information on the aj 
proximations involved in the velocities. 
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Heat Addition to a Flowing Gas 
Roger Weatherston 
Cornell Aeronautical Laboratory, Inc., Buffalo, N.Y. 
November 28, 1949 
HE FOLLOWING RESULTS of experiments that have beet 


conducted at Cornell Aeronautical Laboratory may be ol 
interest in connection with recent discussions'~® of the problem 
of heat addition to a flowing gas 

A 6-in. diameter, 6-ft. long tube was placed axially in a free 
jet of air. Heat was added to the internal flow by means of 
Picatinny Arsenal 13-B flares that discharged transversally into 
the flow through flush wall openings at a short distance (about 
18 in.) from the tube inlet. Inlet and outlet Mach Numbers 
and stagnation temperatures were measured by means of stati 


pressure and stagnation and temperature probes 


The flares provided a steady source of heat for approximately 


pressure 


45 sec., which was ample time to establish a steady régime of 
heat addition and to record instrument readings. 

The experimental results are plotted in Fig. 1, together with 
the theoretical curves calculated by the procedure of reference | 
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with consideration of wall-friction forces. These forces, measured 
in cold flow, were introduced as an additional impulse term in 
8) of reference 1 and were assumed to be constant in the 
computation of the effects of The addi 
tional impulse caused the initial values of 1/4, and M; to be differ 


ent 


Eq 


heat addition. 


The system tested was that of Case II of reference 1 (constant 
stagnation temperature and pressure upstream and constant 
static pressure downstream of the heating region). 

The experimental results were in good agreement with the 
theoretical prediction and proved that important modifications 
of the upstream flow conditions are produced by changes in 
heating rate even when choking does not occur anywhere in 
the flow 

This work was sponsored by the U.S. Office of Naval Re 
search 
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Potential Flow Around a Rotating Cylindrical 


W 


Blade 


. R. Sears 


Cornell University, Ithaca, N.Y 


De 


cember 2, 1949 


THE PROBLEM of an infinitely long cylinder, of 


Consise 
arbitrary section, rotating steadily about an axis normal 


to 
mm 
cy 
ve 
de 


wi 


it, in a perfect incompressible fluid otherwise at rest. This 
iy be thought of as an idealization of the case of a rotor of 
lindrical or nearly cylindrical blades. Let g denote the velocity 
ctor; @, the angular velocity; and p and p, the pressure and 
In addition, since it is convenient to use axes rotating 


nsity. 
s, Wherein y is 


th the blade, let us define a Cartesian set x, y, 


measured outward along the blade from the center of rotation, 
> is measured along the axis of rotation, and x forms the third 


member of a right-hand system 
moves in the 


r 
ve 


an 


pa 


i.ec., for y > O the cylinder 
—x direction. We also define the radius vector 
(x, y, 0), the unit vector k in the z direction, and the relative 
locity q’ gq-—-® Xr. 

We have, then, 


curl g’ = —2wk, div g’ = 0 (1) 
d by integration of the Euler equations for this case! 
(p/p) + (q’?/2) = (1/2)w*r? + constant (2) 


As usual, the solutions of Eqs. (1) are found by combining a 
rticular integral with solutions of the homogeneous problem. 


We write 
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u'’ = dr typ, | rr + buy + be = O 
yy = by — vz where ) Yer + byw = 2w (3) 
w’ = gz ) 


That is, we assume that the particular integral will be expressible 
by means of a two-dimensional stream function, and we propose 
to find the irrotational part of the solution in the form of a three- 
dimensional velocity potential. Here u’, v’, w’ denote the com- 
ponents of q’, and it is easily verified that Eqs. (1) are satis- 
fied. 

The boundary conditions are: 

(1) The normal component of q’ vanishes at the cylinder 


i.e., 


lu’ + nw’ = O on C (4) 


where C is the contour of the cylinder cross section and (/, m, n) 


is the normal to C. 
(2) Asx— mo and/orz— ~,q’— —® X r—ie., 


u’—>wy, v’—>—wx, w’' >0 (5) 


A particular solution satisfying the differential equation for 
y is 

y = (1/2)wr? (6) 
If this value is selected, the boundary conditions become 

(1) l(g@r + wy) + nd: = O on C (7) 


(2) dz, by, dz > Oasx — om and/or s—> «= (8) 


These are exactly the boundary conditions for the plane flow 
produced by motion of the infinite cylinder in the —-x direction 
with speed wy, and the differential equation for the potential of 
that flow is the same as for @¢ if ¢,, = 0. In other words, if 
flow 


¢i(x, z) denotes the potential for plane, steady past 
the cylinder placed in a parallel stream of unit speed, we can 
write 

@ = wylgi(x, z) — x] (9) 


Consequently, combining our results, we now have 


ul = dr + wy = wyei,(x, 2) 
v = dy — wx = wl[¢gi(x, 2) — 2x] (10) 
w’ = $2: = wy¢i,(X, 2) 


This solution satisfies all the requirements of the problem. 

Thus, the flow in a plane transverse to the cylinder at any y 
is the same as steady plane potential flow about the same cylinder, 
so far as the uw’ and w’ components are concerned. The outward 
component v’ is given by a simple formula, dependent upon the 
potential of such plane flow, and is independent of y. 


It may be of interest to look at the special case of the circular 


cylinder. 
Circular cylinder of radius 1: yg, = (s + s~') cos @ where 86 = 
tan-1(s/x) and s = (x? + 2?)' Here, 
v’ = w(s + s~! — 2s) cos 6 = w(s~! — 5) cos 8 
Thus, on the circular cylinder itself, v’ = 0. The flow imme 


diately adjacent to the cylinder lies in a transverse plane. 

It should be noted that we have implicitly assumed here that 
there is no circulation about the cylinder. If the circulation 
were not zero or at least constant along the length of the cylinder, 
the presence of trailing vortices would have to be considered 
The results given here are being applied by L. E. Fogarty and the 
present author to calculate laminar boundary-layer flow on rotat 
ing blades. Subsequently, attention will be directed toward 
the problem of the blade with circulation. 
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f'(xe) =-- e~* sinh (1 — e~®) cosh 6 do = 
x 0 
] 1 e 
_~ P x sinh (9 ee e_ e 20) 469 
x 2 
and since! 
x sinh @ _ I Ni(x fee 
do = ~ [No(x) + Eo(x)] 
0 ~ 
2 ; k 
e7* sinh 6 + °d@ = ’ f 
0 
= )k + 1 
- [Sx(x) = Ne(x) = Ex(x)], 8 > O 
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A Note on the Evaluation of a Definite Integra) 


Henry E. Fettis 


Mathematician, Dynamics Branch, Aircraft Laboratory, Air Mater 
Command, Wright-Patterson Air Force, Dayton, Ohio , 


November 23, 1949 


N THE THEORY of the oscillating wing of finite span, it js; 
quired to evaluate the following integral: 


« 


-tt / 2 2 
, (3 .= 2° ing 
é + Vv + F) dt 
0 st 
e~ #2 sinh 6/4 + sinh 6 — cosh @)cosh 6 dé 


F(iz) = 


0 sinh 6 


Series expressions for this function are given by Cicala! and }y 
Dingle and Kiissner.?_ In the following development, the functig 
is expressed in terms of known functions, with the exception ¢ 
one part, which, however, may be expressed by a simple an 
rapidly converging power series. 

It is convenient to consider first the function with the argument 
real—namely, 
“© -—x sinh (1 4 sinh @ — cosh 8)cosh 6 dé 

F(x) = . 9 
0 sinh @ 

It is known that F(x) has a logarithmic singularity at x = 
and is regular elsewhere and, further, that 


lim [F(x) + log x] = 1 — y — log 2 


x—-0O0 


where y = 0.57732 is Euler’s constant. It is therefore sufficient 
to consider the function 


- 
e~* sinh 677 4 sinh @ — cosh @)cosh 6 dé 
f(x) = log x + 


0 sinh 6 


which is regular at x = 0 and may therefore be developed in 
Taylor Series about this point 
Differentiation of Eq. (4) with respect to x gives 


« 


in which Nx is the Neumann function, Ex is the Lommel-Webet 

function, and S; (k > 0) isa polynomial in (1/x) determining the 

singularity of N;, it is found that 

e l renee ; ; - 

f"(x) = _ 12.51(x) + So(x) + mw [Eo(x) + F2(x)] 4+ 
x 


w[No(x) + N2(x)]} 


Making use of the following known relations: 


2 4 2 { 
Si(x) = -, S(x) =—, Eo(x) + E(x) = E\(x) — 9 
x x? 4 xT 
No(x) + Ne(x) = (2/1) Ni(x) 


gives, finally, 


Therefo 


and sin 


Thus, 


except 


and tl 
pansio 


expan 


In 


pressi 


Fr(iz 


Fy (12 


in wh 


Trai 
is re 


Dis: 
Apr 


bric 








Integra 


Air Mate, 
hio — 


an, it isp 


h 6 d6 


la! and by 
he functig, 
Ception of 
imple an 


argument 


6 db 


itr = 


sufficient 


sh 6 dé 


ed ina 


6 (5 


Veber 


g the 


Be 1 vr : ( 
"x)= ¢ -- ~ E\(x) + Ni(x)] > (9) 
J (x ¥ x? a, | 1 (Xx I 
Therefore 
+ : 
) 1 us E\(x) { 
= f(0) + = <x — 
t) = JM) 0 ix 2 C4 f 
- Os x [ Ni(x) ]\ 
a { dx 
. ” «& a. > (10) 
mi a | E,(x) | te 
= fi(( : — > dx 
wre. ish s | 
E Ni (x) + / = {" x / 
Ni (x) - No(x) dx | 
2 x 2 Jo ) 
and since 
na wx 
No(x)dx = - [ No(x) Ei(x) — Ni(x)Eo(x) ] (11) 
0 as 


wx 
Fix) = 1 — y — log 2x — ; { No(x) FE, (x) — Ni(x) Eo(x)] + 


e 1 “a ] w E\(x) 
Ni(x) + + a . dx (12) 
2 x e Ls 2 & 


Thus, F(x) is determined in terms of known functions with the 


*T1 «BE 2] 
— - | dx 
0 x 2 x 


and this term presents no great difficulty since 
pansion can be written down at once from the corresponding 


exception of the term 
u(x) = (13) 


its Taylor ex- 


expansion of F,: 


Tr E 1 x? 4 x! x6 ; 
(x) = —_ = ; 
hS FS*6 8 86IS*O7 (14) 


In the practical application, x = iz, and the following ex- 


pressions result: 


Feliz) = 1 — y — log 2: + 
Wz T , 
4 [Jo(z)Li(s) — T(z) Lo(s)] — 9 A102) + u(iz) 
‘ T T : . 1 
Fi(iz) = — = + = 2{Ko(z)Li(s) + Kyi(z)Lo(s)] + | Ku = ‘| 


in which 


These expansions are convergent for all values of z and are 
particularly useful when z is small. 
A short table of the functions Lo 
by Owen.5 


and L,; has been compiled 
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Lift-Curve Slope for Swept and Unswept Wings 


Innes Bouton 
McDonnell Aircraft Corporation, St. Louis, Mo 
October 18, 1949 


A’ EXPRESSION has been developed by the Loads Group of 
McDonnell Aircraft Corporation for the low-speed lift- 
curve slope of swept and unswept wings. This expression seems 
to correlate with experimental results over a wider range of 
aspect ratios and sweep angles than any other known formula. 
This expression is as follows: 


. C la 
Cla = : 
1 1 18.24Cie (1) 
7 + , 
V cos A A A 
where Cra is the surface lift-curve slope per degree, Ci, is the 


two-dimensional lift-curve slope per degree of the section per- 
pendicular to the sweep axis, A is the angle of sweep at the 
quarter chord line, and A is the aspect ratio. 

The idea for Eq. (1) resulted from a similar expression de- 
veloped from the Jones edge-velocity correction by George, of 
the Bureau of Aeronautics, who presented the expression: 

Crm 
1 18.24Cie 

4 
A A 


+ (2) 
for unswept wings. The first two terms of the denominator 
represent the ratio of semiperimeter to span for an elliptic wing 
and are empirically applied to a wing of any plan form. The 
same reasoning for a swept wing results in the term 1/cos A + 
1/A. The formula incorporating this term was compared with 
experimental data, and the correlation was found to be poor for 
swept surfaces. An empirical change of the term 1/cos A to 

V cos A, as is often done in dealing with the effects of sweep, 
results in Eq. (1), which correlates satisfactorily with over 90 
surfaces ranging in aspect ratio from less than 1 to over 12 and 
with more than half of the surfaces having sweepback, some as 
much as 60°. (1) and 
experimental data is no more than the experimental differences 
A value 


The maximum departure between Eq 


between similar surfaces or approximately 10 per cent. 
for Cig of 0.105 is suggested when the section lift-curve slope is 


not readily available. 


Base Pressures at Supersonic Velocities 


Freeman K. Hill 
Physicist, Applied Physics Laboratory, The Johns Hopkins University, 
Silver Spring, Md 
November 21, 1949 
N AN EARLIER ISSUE of the JOURNAL OF AERONAUTICAI 
ScIENCES,! data were presented of pressure measured at the 
base of missiles in free flight in the range of Mach Numbers 
from 0.9 to 2.5. The reduced data in terms of base pressure 
coefficient indicated substantial differences with data obtained 
by other investigators and by other experimental techniques. 
However, since the time of this publication, the free-flight data 


have been substantiated by extensive wind-tunnel tests with 
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Fic. 1. Base pressure coefficient versus Mach Number. The APL/JHU free-flight data and all the wind-tunnel data shown p. 63 
were taken on a model having an ogival nose tangent to a cylindrical body, four biconvex tail fins, and a flat base. The length-to. Th 
diameter ratio was 6.1. The OAL wind-tunnel data were corrected where necessary to compensate for the effect of reflected shocks 


from the tunnel walls. page 

homologous models at relatively large Reynolds Numbers, however, are so predominant that the base pressure coefficient 

approaching but not equaling those of free flight. In addition, cannot be simply expressed in analytical form. 

a recent test at a Mach Number of 5.1 has provided data that For relatively high Mach Numbers we have therefore con 

enable a correlation to be made of nearly all results at Mach sidered® that a fraction, f, of the base pressure, p», is proportional 

Numbers above 2 and also suggest the causes of the wide dis- to the dynamic head, (1/2)p,V*. This yields the dynamical 

crepancy that exists in many of the base pressure measurements equation 

at low supersonic Mach Numbers. It should also be mentioned . bi 

that the recent letter of Gabeaud to the JouRNAL? further sub- dp/p, = (1/2)AV? | On 

stantiates the general analysis presented below. (1/f) pb Li 
Measurements were made at M = 5.1 in the 5- by 8-in. test When this equation is integrated and the results are written it 

section of the APL/JHU blow-down wind tunnel. The test terms of the base pressure coefficient Cz, one obtains Arth 

model was a 1/3.2-scale model of one of the free-flight models , eo P Asso 

described earlier,! and the test Reynolds Numbers were 5.4 to Se m po — Br = : ( Poe eee aur)’ ; nee. r Ur 

8.5 millions. The model was supported in the wind tunnel by q yM? 2 yi vece 


extending two of the tail fins to the test section walls, and the é 


pressure was measured at the center of the base. Numerically where g; = (1/2)p,V2, and M is the local Mach Number. Eq 


the base pressure coefficient, Cz, was found to be —0.0457 with (2) now presents two arbitrary constants, \ and f, which may be | 
a standard deviation of 0.0009. (Complete details given in gyalyated from the experimental data. The experimental Dist 
reference 3.) Cg is defined by Ca = (pp — pi)/qi, where po is values of Cz found at Wo = 2.50 and 5.10 were —0.132 and M. | 
the absolute pressure at the base of the missile, p; is the free- —0.0457, respectively, and substituting these in Eq. (2) gives wae 
stream static pressure, and gq, is the dynamic head A = 0.055 and f = 0.5394. Eq. (2) is plotted in Fig. 1. Not mat! 
Fig. 1 is a summary of base pressure coefficient data now only does this curve show the trend of base pressures at high has 

available from both free-flight and wind-tunnel tests made with — ygach Numbers, but also in the Mach Number range near 2 its subj 
homologous models having tail fins. The free-flight data are steep descent indicates where the effects of Reynolds Number in di 
reported in reference 1, and the wind-tunnel data were obtained boundary layer, body shape, etc., play a predominant role. Fo for ¢ 
in the APL/JHU tunnel at M = 5.1 and in the Ordnance Aero example, at M = 1.75 the knee of the free-flight base pressurt _ 
physics Laboratory wind tunnel, Daingerfield, Tex., at 7 = 1.5, curve indicates the region where divergence from the dynamica whic 
1.75, and 2.00. Free-flight results obtained by ballistic firings! PETE OG be e 
are also shown, where the model configurations and test condi It may be noted from Eq. (2) that for values of M > 95, ~ 
tions were different from those mentioned above; these results, p» becomes negative, and, hence, Eq. (2) is no longer physically the } 
however, illustrate the dispersion of data in the low supersonic tenable. However, at these high Mach Numbers the bast brief 
range of Mach Numbers and the trend toward agreement above pressure is going to be nearly zero, and there seems little poial Bi 
M = 2. The theoretical curve of Gabeaud is also shown. The in deriving an analytical expression for it, especially in view of ve 
agreement in the data above JJ = 2 is evident from this figure. the fact that the usual isentropic laws are questionable in flows who 
In this connection it may be pointed out that as V/ increases, where such high Mach Numbers obtain. For the purpose of dire 
the value of the base pressure will approach but cannot exceed predicting base pressure on simple missiles in the range 0.9 « is tl 
the perfect vacuum condition, which, according to the above M < 1.8, experimental data must be used, but in the range - 
definition of Cs, g ves a value of —2/yM?. It is clear, therefore, 2 < M < 6 an expression such as Eq. (2) is useful. 
that at high Mach Numbers the base pressure becomes mainly 
a dynamical effect and will approach the perfect vacuum limit. REFERENCES — 
This suggests that an expression may be derived similar to the 1 Hill and Alpher, Base Pressures at Supersonic Velocities, Journal o 
one set forth by von Karman and Moore in 1932 and that the the Aeronautical Sciences, Vol. 16, No. 3, p. 1538, March, 1949 ; 
experimental data now available may be used to fix the boundary Tee, See oN sips 624 ag nines hy ecpeey Pcagencaanl ~~ . 

: 2 Journal of the Aeronautical Sciences, Vol. 16, No. 10, p. 638, October, 1949 von 
conditions. In the low Mach Number range, 1 < 7 < 2, the Hill. Bose Pressures of a Mach Number of 5:1, APL/JHU Intern the 


effects of model shape, Reynolds Number, boundary layer, etc., | Memorandum CF-1306, July, 1949 
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Erratum 
William F. Savage 
Assistant Professor of Mechanical Engineering, University of 


Kentucky, Lexington, Ky. 
December 30, 1949 


T® FOLLOWING ERROR appeared in my Readers’ Forum item 
“Effect of Disposable Load on Range of Turbojet Airplanes”’ 
JoURNAL OF THE AERONAUTICAL SCIENCES, Vol. 16, No. 10, 


p. 637, October, 1949). 
The first equation at the top of the right-hand side of the 


page should be: 


: [Kr/(1 + K')] 
Vave = 4V; = > , 9 
;' (Zs + OK + ey") 
"LK\ 2K’? +2 — Kr 
- + 


On the Surface Effects of a Compressible 
Laminar Boundary Layer 


Arthur N. Tifford 

Associate Professor of Aeronautical Engineering, The Ohio State 
University, Columbus, Ohio 

December 5, 1949 


HE PAPER, ‘“‘Temperature and Velocity Profiles in the 
Compressible Laminar Boundary Layer with Arbitrary 
Distribution of Surface Temperature,’ by D. R. Chapman and 
M. W. Rubesin (published in the September, 1949, issue of the 
SCIENCES) an excellent 
A notable contribution 


JOURNAL THE AERONAUTICAL is 
mathematical analysis of the problem. 
has been made to the analysis of the laminar boundary layer 


I find myself 


OF 


subjected to an imposed temperature gradient. 
in disagreement, however, with the emphasis given to the need 
for considering the change in the velocity profile due to tempera- 
ture effects when dealing with most of the practical problems to 
which reference is made in the subject paper. Rather, it is to 
be expected on the basis of earlier works that the Blasius velocity 
profile persists effectively to relatively high Mach Numbers 7/ 
the fluid properties are taken as those at the plate’s surface. A 
brief review of the relevant literature follows: 

Busemann in his original treatment! of the laminar boundary 
layer along an insulated flat plate exposed to the flow of a gas, 
whose viscosity and thermal conductivity vary with temperature 
V R, increases only 20 per cent 


ut Car R, 


directly as VY 7, noted that « 
isthe Mach Number increases from 0 to is expressed 


as 


Mos (Ou oy (pu ll 2/2 Vv U 


rather than ; 


is 


My (Ou Oy) (Pa U «~? 2) OFF Vo 
von Karman, in his famous Volta Congress paper,” suggested 


the application of this finding to the analysis of the effect of 
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compressibility on the turbulent boundary layer. However, 
in the later paper* by von Karman and Tsien reanalyzing Buse- 
mann’s problem, assuming the coefficient of absolute viscosity 
and tirmal conductivity of the gas to vary directly as 7°", 
no consideration was given to this remarkable simplification and 
its implications. Thus for Eckert and Drewitz* 


to show that von Karman and Tsien’s results confirm Busemann’s 


it remained 


finding of only a slight dependence of ca; R, (properly defined) 
on the Mach Number. In addition, they find that only a small 
deviation from the Blasius velocity Mach 
Numbers as high as those considered in the paper being discussed 
here—i.e., Mach Numbers up to 5—when the nondimensional 
boundary-layer coordinate is based on the kinematic viscosity of the 
gas at the surface of the plate. What are the independent findings 
of the present paper? 
Eq. (58a), 


profile occurs at 


1.328 
= 7 vV " 

V R, 
was derived as applicable to a plate at a constant temperature 
Direct substitution of the definition of C and R, and simplifica- 


Tw 
Cay = aged 


poi’ 


tion yields 


(1) 


(te V Ucl./ve)/ (ow Ua?/2) = 1.328 
exactly the result obtained by ignoring velocity profile and fluid 
property alterations due to the temperature field! Furthermore, 
it is the comparison in Fig. 6 between Eq. (58a) and Crocco’s® 
more exact results that was used to establish the precision of the 
new analytical approach as within 5 or 6 per cent for Mach 
Numbers as high as 5. This identical precision therefore applies 
to an analytical approach neglecting fluid property variations. 
Fig. I compares lines representing Chapman and Rubesin’s 
CT*, with plots of the actual variation of 
It is seen that 


* = 


approximation, yu 
the absolute viscosity of air with temperature. 
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only at extremely low free-stream temperatures is the approxi- 
mation reasonably good. Under these conditions, Cy is essentially 
the same as C. Eq. (58) for the total drag coefficient, ca,, 


f"0) | 7! Codx* 


“LV RC 0 Vx" 


therefore reduces for extremely low free-stream temperatures to 
Eq. (I) above. Jn other words, when a plate has a nonuniform 
surface temperature, an analysis based on constant fluid properties, 
properly interpreted, yields the same results as the analysis under 
discussion, whenever the latter’s assumption, u* = CT™, ts actually 
a good approximation. The poorer the approximation becomes, 
the less the agreement between the two analyses. Which is the 
better approach can only be fully established after a more exact 
analysis of the boundary layer of a plate of nonuniform tempera- 
ture has been made for atmospheric flight conditions. 


(58) 


The discussion up to this point has been primarily concerned 
with the surface friction. One may readily confirm, however, 
by an analogous study that the heat transfer rate determined by 
an analysis neglecting fluid property variations is identical with 
that derived in the paper under discussion, except when dealing 
with a plate of nonuniform temperature operating in a flow 
field having a free-stream static temperature of above approxi- 
mately —100°F. 
cussed above, indeterminate at the present time which type of 
analysis is the more precise. 


Under the latter circumstance it is, as dis- 


It may be of interest to the reader to replot the boundary-layer 
velocity profiles of the typical example of the paper as a function 
of (y/2)V Uw /vyx instead of (y, av ‘o/Vex. The three 
distinct velocity profiles shown in Fig. 10 merge into a single one, 
the Blasius profile, except for an extremely small region near 
the outer edge of the boundary layer. This is not surprising in 
view of the discussion which has been given here. 


In conclusion, (a) Chapman and Rubesin’s paper is a major 
contribution in that it time the 
solution of the boundary-layer equations in the presence of a 
fairly general analytic distribution of surface temperature. 


makes available for the first 


(b) The results of the analysis, properly interpreted, yield 
additional support for the assumption that the evaluation of local 
botindary-layer surface phenomena, such as friction and heat 
transfer in a compressible gas flow field, need not be complicated 
by the consideration of the fluid properties anywhere but at the 
surface if the Mach Number outside the boundary layer is less 
than about 5. 
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On the Vorticity of the Flow Behind Curve; 
Shock 


S. |. Pai 


Institute for Fluid Dynamics and Applied Mathematics, Univers 
of Maryland, College Park, Md. ; 


December 1, 1949 


[ IS WELL KNOWN that vortices are generated suddenly }y 

shock waves and that the strength of the vorticity is pes 
portional to the rate of change of entropy with respect to streap 
function.'. One of the assumptions usually made in calculatio, 
of flow involving shock wave is that the vorticity introduc 
into the flow by the shock wave is extremely small and ne 
ligible. This note will show that this assumption, in gener,j 
is not accurate and that, in most practical cases, the vorticity 
introduced by the shock should be considered, particularly whe 
the flow behind detached shock is considered. 

For simplicity, we consider only the practically importay 
iso-energetic flow. The basic equations for the present analysi 
are :? 


p .ds 
R’ & 


So —- S 3 inf 
= log ( : M,? sin? a — : ) 
C, y¥+ 1 y+ 1 


y¥-1 2 Y 
a 
y¥+1 (y + 1)M;,? sin? a 


where p is the pressure, FX is gas constant, y is radial distance, 


vorticity = w = 


e = 0 for two-dimensional flow, « = 1 for axially symmetrical 
flow, S is entropy, y is stream function, C, is specific heat at 
constant volume, y is ratio of specific heats (C,/C,), M is Mach 
Number, a is shock angle with respect to velocity in front of 
shock, sub 1 refers to the value in front of the shock, and sub? 
refers to the value behind the shock. 

In order to show the order of magnitude of vorticity introduced 
by the shock, we consider the small perturbation from those 
irrotational flow involved shock wave, such as a uniform super- 
sonic flow over two-dimensional wedge or over unyawed cone 
There are two possible kinds of perturbation from the original 
irrotational flow: 

(A) Disturbance comes from the front of the shock, which 
affects both MM, and a. 
case of such a problem.* 

(B) 
affects only a. 
the body is slightly different from that of the wedge or cone 


For instance, Adams studied a specia 


Disturbance comes from the rear of the shock, which 
This corresponds to the case that the shape o 


Since case (A) involves both the change of J/, and a, we need 
to consider it only. The change of entropy behind the shock 


with respect to stream function due to the small perturbation is 


dS» 1 OS: Oa OS; OM; 
-|. : 
re) A 1; Oo 


Oa Oa 
where p is the density, g is the resultant velocity, 6 is the shock 


dy pogoy® sin 6 


angle with respect to go, and o is arc length along the shock curv 
S; = constant. 

The first order value of vorticity behind the shock introduced 
by the perturbation is 


‘ p» QS2 da’  ‘dS2 OM;’ 
o= - mu 3 , 4 
pogeR sin d| 0a Oo 0M, Oo 
where prime refers to the value of perturbation quantity. 
1 OS: _ 2yM;? sin 2a[(y — 1) M,? sin? a + 2] — 47 cot a[2yM,? sin? a — (y — 1)] 
Cy 0a [(y — 1)M)?2 sin? a + 2][2yM)? sin? a — (y — 1)] 
and 
1 OS: a 4y(y — 1)(M? sin? a — 1)? p 
C, OM, Mil(y — 1)M,? sin? a + 2][2yM,2 sin? a — (y — 1)] 
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The ne pe curves of (1/C,)(OS2/Oa@) and (1/C,)(0S2/0M;) 
at M, = 3 against the shock angle @ are plotted in Fig. 1, which 
shows that these quantities are of the order of magnitude of unity 
over practically the entire range of shock strength. The first 
term of Eq. (4)—i.e., (0S:/0a)(Oa’/O00)—shows the variation of 
vorticity along a given shock polar (MM, = constant), and it is 
small only near a = Mach angle and a = 90°. The second 
term of Eq. (4)—i.e., (0S2/0M,)(0M,'/0a)—shows the variation 
of vorticity from one shock polar to another, which is only 
negligible when @ is near Mach angle. 

In general, the vorticity introduced by the shock wave is of 
the same order of magnitude as the perturbation terms 0a’/d0¢ 
and 0M,’/dc, and it cannot be neglected. In case of detached 
shock, the shock angle varies from normal shock a = 90° to 
Mach angle a = sin~'(1/M,); it is clear that the vorticity 
introduced by the shock should be considered there. 
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Formulation of the Flutter Problem for 
Solution on an Electronic Analog Computer 


Eugene F. Baird and Henry J. Kelley 
Grumman Aircratt Engineering Corporation, Bethpage, L.!., N.Y 
November 93, 1949 


HE PRESENT NOTE deals with the formulation of an initial 
value problem in the stability of airfoil oscillations and its 
solution by means of an electronic analog computer. The 
advantages of this approach lie mainly in speed and ease of 
analysis 
Wagner!’ first obtained the lift response as a function of time 
for an airfoil whose quasi-steady circulation I" experiences a 
sudden unit increment. In terms of lif€ deficiency ¢: 


I J (V(t i - i (t d l 
.= J “t) — -_ ) (1) 
' { Jo a " b . "| 


Using the notation of Theodorsen,? the circulatory lift may be 


Written: 





FORUM 189 


t f V 
L = 2xpVb 41 - o[ (t — nt x 
0 | b 


: 
lh+ Va+t+b('/2 — a)a|dr (2) 
T 


That the quasi-steady circulation depends on the downwash at 
the rear quarter chord in the above fashion has also been shown 
by Kiissner.’ 

The equilibrium of forces and moments is expressed in the two 
following equations, which include virtual inertia terms: 


Mh + Sgt + Ch + mpb*(i + Va — aba) + 
t , 
V \d 
2brpV J, _ (t — | x 
an of | { dt 
[h + Va + b('/2 — ala] dr = 0 (3) 


Tqa + Sah + Caa + arpb?{('/2 — a) Vba + b7('/s +a%)a — abh| = 


V { 
2rp Vb2("/s a) {1-6 (t — | x 
0 b f 


4 E + Va + b('/2 — aya fer = (0 (4) 


Since the function @ is not available in closed form, various 
approximate representations have been developed. A number 
of these may be listed: 


o(c) ~ 0. 165¢e—9 04550 4 0.335e—9 300 (5a)* 
o(o) & 0.25e— 9-5" 4+ 0.25(1 + 0.1850)e—9 189" = (5b) 
o(o) & 0.165¢e—9-% 1 + 0.335¢e—9-32¢ (5c)* 


For comparison with Theodorsen’s analytic result, the Fourier 
transforms of 1 — @ in the approximate forms above have 
been plotted in Fig. 1. The scale is chosen for purpose of con- 
trast. 

Choosing the Jones approximation for ¢, we rearrange as 


follows: 


o(o) & ('/2) — [0.165(1 — e— 994557) + 0.335(1 — e— 9-302) | 


The method for treating the integrals in Eqs. (3) and (4) is 
to introduce auxiliary variables — and n into the equations of 
motions, creating a set of four simultaneous linear differential 


equations in four dependent variables as follows: 


(M + wpb2)h + (Sa — mpbia)a + Crh + wpb? V(3/2 — ala + 
mpbVh + rpbV2a +E +n = 0 (6) 
[Ta + mpb*('/s + a?) la + (Sa — mpb'a)h + Caa + 
mpb3(!/, — a)?Va — pb? Vso + adh — rpb?('/2 + a) X 
V2a — d/o + all—E+n) = 0 (7 


b ‘ 
+ ). = (2rpVb)(0.165)[h + Va + b('/2 — aja] 
0.0455 V 


(8) 
b F * s : 
n + . 7 = (2rp Vb) (0.335) [h + Va . bt} e@=~ @2ja@ (9 
0.301 


It is easily shown by Laplace transformation that £ and » 


represent the integrals: 


7 


£ = (2xpVb)(0.165 / [1 — ¢—0-:0655(V/0)E — 1) x 
/0 
a4 : rca —— 
(fh + Va + b(!/2 — alaldr + (2xpVb)(0.165) X 
(1 — e~ 0.0455 Vit/b)) h + Va + b('/e — aaj 0 
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Solution of linear differential equations with constant 
coefficients is a routine task for most electronic analog computers. 
The aerodynamic forces and moments that arise from initial 
values of h, h, a, & are provided for in the representations £ and 
n. The procedure consists of solving the four equations at 
successive values of forward velocity for sets of initial conditions 
chosen with the purpose of exciting pure modes in each of the 
degrees of freedom. For values of V greater than a flutter 
however, arbitrary initial conditions suffice to 


In practice it has been found that use of V?/b as inde- 


speed, excite 
flutter. 
pendent time variable facilitates the analysis in that only a few 
coefficients depend on V. 

Results of a two-dimensional two-degree numerical example 
will be shown for comparison with results of a conventional 
(A.M.C. Method). The following parameters 
0.969 ft.; a = 0; 
= 305 rad. per sec.; 


were 
natural frequencies w, = 147 
rad. per sec. and wae M = 15 lbs. per ft.; 
Sa = 0; Iq = 675.7 Ibs. ft.? per ft.; p at 10,000 ft. altitude. A 
plot of damping coefficient g as obtained by the A.M.C. method 
appears in Fig. 2, together with rates of decay observed on the 
R.E.A.C.* in the form of logarithmic decrements. Although 
flutter speeds are identical, the steady state assumptions of the 
A.M.C. method preclude comparison for other V.7 

The method affords facility in the study of gust response and 
in the inclusion of the effects of servo units in control surface 


analysis 
used: b = 


analysis. Extension to the supersonic case does not follow di- 
rectly, however, since the response functions involved are not 
easily represented as solutions of linear differential equations 


(discontinuous derivatives). 


* Reeves Electronic Analogue Computer 
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Boundaries of Fluid Mechanics 


Research Associate, Aeronautical Research Center, 
Michigan, Ann Arbor, Mich 
November 14, 1949 


Roberts points out < 
Liepmann to distinguish the realms of fluid mechanics. In the 
attempt to define (4) 


BETWEEN CONTINUUM ANALYSIS 


AND SLIP FLow 
boundary layer thickness 


Reynolds Number 


viscosity coefficient 
= mean molecular velocity for air 
= velocity of sound 


Subscript 0 signifies sea-level conditions. 
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Comparison of slip flow criteria. 


(2) Tsien: 


»/s = M/VR, R>1 
\/s = M/R, R<1 


University of 
4/6 = O.O1 


(3) Schultz, Spencer, and Reifman:! 


Review, Mr. 1/6 = 0.4388(M/V R) = M 2VR 
, criterion! suggested by H. W. »/6 = 0.05 


Roberts:! 


compares these R/M2 = 100 


We can compare these criteria in two ways: first, 


R plv/p 7 l pa 7 l Ao 


In the region of the atmosphere under consideration 7/7) can be 


considered unity in an order of magnitude analysis. 


M2 Mar (7 i. My 1 
= (0.4388) = 
R l To 1 2 


Second, an extension of Donaldson’s Appendix A, 


M? M? M? 
R pul /p (yl d) (pr/n) 
but p» = 0.499p cd and c = 1.462a; hence, w« = 0.73par and 
M? M? " , Mr 3 Mx 
= ().7 ~ 


R_11/0.73.a . “4 





M M (v/a) M uw 7 Md po( T/T )9-78-0-5 
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Let us use M2/R = (MX/I)(1 2) to compare the criteria. 
(1) Donaldson: 
MyX/l = 0.04 
(2) Tsien: 
MyX/l = 0.0002 
(3) Schultz, Spencer, and Reifman: 


My/i = 0.02 


(4) Roberts: 
Mx/l = 0.02 


It thus appears that in any order of magnitude discussion that 
criteria (1), (3), and (4) are equivalent. 

Lo® has superimposed Tsien’s criteria on a graph of WV versus 
log R for the flight of V-2 missile No. 21, fired on March 7, 1947. 
We have added the other criteria to this curve. We note in 
Fig. 1 that there is a spread of 131,000 ft. between application of 


these criteria. 


FOR THE LOWER BOUND OF FREE 


MOLECULE FLow 


DIFFERENT CRITERIA 
(A) Roberts: 
R/M? = 0.01 


M/V/R = 10 
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(B) Schultz, Spencer, and Reifman: 


\/6 = 10 
\/5 = 0.438(M/V R) 


Tsien: 
\/6 = d/l = M/R = 10 

One notes that, in the intended region of application of criterion 
(B), it was thought that R > 1, while criterion (C) implied 
that R< 1. Thus we can make no real comparison betweeg 
these criteria. 

Of course, we do not know which criteria are the best, sings 
this answer must be determined by experiment. 
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The Pressure on a Slender Body of Nonuniform 
Cross-Sectional Shape in Axial Supersonic Flow 


(Continued from page 175) 


Fig. 3 shows a sketch of the distorted element, and 
Fig. 4 gives the ordinates and pressure coefficients for 
lines on the surface defined by various planes passing 


through the flow axis. The importance of the longi- 
tudinal shape of the element is indicated by the pres- 
sure reversals that occur for 6 = 0° and@ = 90°. The 
pressure reversals also indicate qualitatively that a 
drag force acts on the forward half of the element and 
that a negative drag force acts on the rear half. This 
tendency to eliminate the resultant drag is consistent 
with the conclusion previously stated—that no result- 
ant drag could be predicted, within the limits of accu- 
racy of the theory, for a body ending in a cylinder. 
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